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Abstract 

Modal analysis provides a powerful tool for efficiently 
simulating the behavior of deformable objects. This 
paper shows how manipulation, collision, and other 
constraints may be implemented easily within a modal 
framework. Results are presented for several example 
simulations. These results demonstrate that for many ap- 
plications the errors introduced by linearization are ac- 
ceptable, and that the resulting simulations are fast and 
stable even for complex objects and stiff materials. 
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1 Introduction 

Interactive modeling of deformable objects has a wide 
range of applications from surgical training to video 
games. Many of these applications require realistic, 
real-time simulation for complex objects. Unfortunately, 
the most straightforward simulation methods turn out to 
be prohibitively expensive for modeling objects of even 
modest complexity. When the high cost of simulation 
couples with the reality that CPU cycles must be shared 
among many tasks, the need for faster, more sophisticated 
simulation methods becomes clear. 

Recently several ingenious techniques for modeling 
deformable objects have been proposed. Examples in- 
clude multi-resolution representations that avoid wasting 
time on irrelevant details {e.g. [4, 6, 8]), reformulating the 
dynamics to make them more stable {e.g. [15,20]), exten- 
sive precomputation to minimize runtime costs {e.g. [9, 
10, 19, 22]), robust integration schemes that afford large 
time-steps {e.g. [3]), and many other approaches that we 
cannot list here due to space constraints. As of yet, none 
provides a perfect solution that satisfies the requirements 
for all interactive applications. 

This paper reexamines a technique known as modal 
analysis that was originally introduced to the graph- 
ics community over a decade ago, but has since been 
largely neglected, with only a couple of notable excep- 



Figure 1 : This example demonstrates a complex model 
being deformed using a modal simulation method. The 
object furthest from the viewer shows the undeformed 
configuration. The nearer objects are being deformed by 
a force indicated by the blue arrows. 

tions (e.g. [10,22,23]). Like the techniques mentioned 
above, modal analysis does not provide a perfect solution 
for every interactive application, but it does provide a so- 
lution that suits some applications quite well. 

The results presented here show that modal analysis 
can be used effectively to model situations where the de- 
formable object is directly manipulated using constraints 
and where it interacts with an environment through con- 
tact forces. We demonstrate that although linear modal 
analysis does incur errors because of the inherent lin- 
earization of the dynamics, these errors are acceptable 
in many contexts, particularly when exaggerated cartoon- 
like deformations are desired. While precomputing the 
modal decomposition for a complex object may take up 
to a few hours of precomputation, for applications which 
make use of fixed content this computational cost only 
occurs during content development and it is well worth 
the dramatic increase in runtime performance. 

The concepts required to manipulate the modal equa- 
tions are to a certain extent conceptually difficult to work 
with but their implementation is surprisingly simple. The 
results shown in this paper (e.g. figure 1) were gener- 
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ated using an implementation that we have ported to sev- 
eral platforms: SGI IRIX, Windows, Linux, and Sony 
PlayStation2. On each of these platforms we were able 
to obtain interactive simulation times even for relatively 
complex models. 

2 Background 

Modal analysis is a well established mathematical tech- 
nique that has been used extensively in mechanical, 
aerospace, civil, and other engineering disciplines for 
several decades. To a large extent the work we present 
in this paper follows as direct application of the methods 
developed in those fields to the task of interactively simu- 
lating deformable solids. There are, however, some issues 
that are unique to interactive simulation, such as impos- 
ing manipulation constraints and computing fast collision 
responses. This paper focuses on those issues. A discus- 
sion of modal analysis and its use with the finite element 
method can be found in the text by Cook, Malkus, and 
Plesha [5], and a more detailed discussion of modal anal- 
ysis, its mathematical theory, and its applications may be 
found in the text by Maia and Silva [13]. 

Modal analysis was first introduced to the graphics 
community in 1989 by Pentland and Williams as a fast 
method for approximating deformation [19]. They used 
a hybrid framework, previously described by Terzopou- 
los and Fleischer [24], that separated the motion of a 
deformable solid into a rigid component and a deforma- 
tion component. The deformable component existed in 
a non-inertial reference frame that moved with the rigid 
component. To avoid the cost of computing the modes 
for a particular object Pentland and Williams used linear 
and quadratic deformation fields defined over a rectilin- 
ear volume instead of the object’s actual modes and then 
embedded the object within the region in a fashion simi- 
lar to a free-form deformation. Although using approxi- 
mated modes is computationally inexpensive, it only gen- 
erates reasonable results for compact objects that are well 
approximated by a rectilinear solid. Pentland and his 
colleagues also integrated their modal deformation tech- 
niques into a interactive modeling system [18]. 

In 1997 Stam developed a modal method for model- 
ing trees blowing in the wind [23]. Rather than starting 
with a deformable object, he computed the low-frequency 
modes from an articulated structure that described the 
tree. Once the closed-form solutions for each mode were 
computed, the response of the tree to a stochastic wind 
field could be computed efficiently. 

Most recently, James and Pai implemented a system for 
computing real-time modal deformations on commodity 
graphics hardware [10]. They focused on modeling de- 
formable skin and soft tissues attached to moving charac- 


ters or as background elements in a surgical simulation. 
Shen and his colleagues have demonstrated an interac- 
tive system that could simulate models with over 10,000 
vertices on a laptop PC with no special hardware acceler- 
ation [22]. 

Other related work includes sound generation tech- 
niques that make use of modal synthesis, and deforma- 
tion techniques that use global shape functions that have 
some general similarities to a object’s mode shapes. Van 
den Doel and his colleagues have used both analytically 
computed modes for simple geometric shapes and sam- 
pled modes from real objects to compute realistic sounds 
for simulated environments [26, 27, 28]. O’Brien and 
his colleagues developed similar techniques that used nu- 
merically computed modes from a finite element descrip- 
tion of an object [17]. Examples of deformation tech- 
niques using global shape functions include: free-form 
deformations and their dynamics extensions [7,21], de- 
formable superquadrics [14], and the boundary element 
method [9]. Modal bases have also proven to be an ef- 
ficient way to compactly encode both shapes and defor- 
mations [11, 12]. Finally, this paper focuses primarily 
on integrating manipulation and contact constraints into 
a modal framework, and there is prior work on applying 
these types of constraints to flexible body simulations [2] . 

3 Methods 

The mechanical properties of an object can generally be 
captured by a function that maps the state of the object to 
a distribution of internal forces. For nearly any non-trivial 
system this function will be nonlinear and the represen- 
tation of state will require many variables. Consequently, 
modeling the object’s behavior over time will involve in- 
tegrating a large, nonlinear system of differential equa- 
tions. These systems are typically far too complex to be 
solved analytically, so some type of numerical solution 
method must be employed. 

Modal analysis is the process of taking the nonlinear 
description of a system, finding a good linear approxima- 
tion, and then finding a coordinate system that diagonal- 
izes the linear approximation. This process transforms a 
complicated system of nonlinear equations into a simple 
set of decoupled linear equations that may be individually 
solved analytically. 

The main benefit of this modal approach is that the be- 
havior of the system can be computed much more effi- 
ciently. Because each of the decoupled equations can be 
solved analytically, the stability limitations that plague 
numerical integration methods are eliminated. Further, 
one may examine each of the decoupled components and 
discard those that are irrelevant to the problem at hand. 
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Figure 2: Using a linear formulation to model a bend- 
ing bar produces acceptable results for small to moderate 
amounts of deformation. For larger deformations signif- 
icant amounts of distortion appear. This example shows 
the deformation corresponding to the bar’s second trans- 
verse mode. 

There are also two drawbacks to a modal approach. 
First, linearizing the original nonlinear equations means 
that the solution will only be a first order approximation 
of the true solution. How objectionable the lineariza- 
tion error is depends on the application and the extent 
to which the objects deform from their initial configura- 
tions. As illustrated by figure 2, small to moderate defor- 
mations exhibit little or no noticeable error when casually 
observed. Even when the errors do grow noticeable, they 
have a cartoon-like, exaggerated appearance that may ac- 
tually be desirable for some applications. 

The second drawback arises because decoupling the 
linear system requires computing its eigendecomposition. 
However we do not believe that this drawback is partic- 
ularly significant. The content in most interactive appli- 
cations is constant, so that eigendecompositions can be 
precomputed during content development and stored with 
the objects. Furthermore, the linear systems are sparse, 
so that fast, robust, publicly available codes may be used 
to efficiently compute the decompositions (e.g. TRLAN 
[29]). 

The remainder of this section describes how one com- 
putes the modal decomposition for a given object and 
how that decomposition can be used to efficiently model 
the object’s behavior. Some of this material has been pre- 
sented elsewhere by others in the graphics community 
(e.g. [10, 19]) but we include it here for completeness. 
The discussion will focus in particular on including ma- 
nipulation and collision constraints in the modal frame- 
work. An overview of the entire process is shown in fig- 
ure 3. 

3.1 Modal Decomposition 

The modal decomposition of a physical system begins 
with a linear set of equations that describe the system’s 
behavior. In general, the equations describing the system 



Figure 3: This diagram illustrates both the preprocessing 
steps used to construct the deformable modal model for 
an object, and the processes that subsequently generate 
interactive motion from this description. 

may be nonlinear, and one obtains the linear equations by 
linearizing about some point, typically the rest configu- 
ration of the system. The linearized equations have the 
general form: 

KdFCdFMd = f , (1) 

where K, C, and M are respectively known as the sys- 
tem’s stiffness, damping, and mass matrices, d and f re- 
spectively as the vector of generalized displacements and 
forces, and an overdot indicates differentiation with re- 
spect to time. The physical meaning of the generalized 
force and displacement vectors, and the method for com- 
puting the system matrices will depend on the type of 
method used for modeling the system. For general fi- 
nite element methods, we refer the reader to the excellent 
text by Cook, Malkus, and Plesha [5]. We are using an 
implementation of the piecewise-linear tetrahedral finite 
element method described by O’Brien and Hodgins [16]. 
Details on computing the system matrices appear in [17]. 

Modal decomposition refers to the process of diago- 
nalizing equation (1). The most general form of modal 
decomposition can be used for nearly arbitrary systems, 
but the systems arising from the finite element method 
we use have a structure that makes them amenable to a 
simpler manipulation provided we assume that the damp- 
ing matrix, C, is a linear combination of the K and M. 
This restriction is known as Rayleigh damping, and al- 
though it is a restriction it still produces results superior to 
the simple mass damping that is most commonly used in 
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graphics applications. With these conditions, diagonaliz- 
ing equation (1) becomes equivalent to solving a general- 
ized symmetric eigenproblem with symmetric, positive- 
definite matrices. Cook, Malkus, and Plesha describe the 
process in detail and we only repeat the end result here. 

With the restriction of Rayleigh damping equation (1) 
may be rewritten as: 

K{d + aid) + M{a2d + d) = / , (2) 

where ai and 0^2 are the Rayleigh coefficients. Let the 
columns of W be the solution to the generalized sym- 
metric eigenproblem Kx ^ XMx = 0 and A be the 
diagonal matrix of eigenvalues \ then equation (2) may 
be transformed to: 
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Figure 4: The two rows show a side and top view of 
a bowl along with three of the bowTs first vibrational 
modes. The modes selected for the illustration are the 
first three non-rigid modes with distinct eigenvalues that 
are excited by a transverse impulse to the bowl’s rim. 


A{z + aiz) + {a 2 Z + z) = g , (3) 

where a = is the vector of modal coordinates 

and g = W'^ f is the external force vector in the modal 
coordinate system. 

Each row of equation (3) corresponds to a single scalar 
second-order differential equation: 

XiZi + {aiXi + a2)zi Zi = gi . (4) 

The analytical solutions to each equation are 

Zi = -h C2e^^i (5) 


where ci and C 2 are arbitrary (complex) constants, and uji 
is the complex frequency given by 


^ — (ofiAi + 0^2) ± \/ {c^iXi + 0^2)^ ~ 4 Xi 

^ • ( 6 ) 


The absolute value of the imaginary part of uJi is the fre- 
quency (in radians/second, not Hertz) of the mode, and 
the real part is the mode’s decay rate. In the special case 
where the term under the radical in equation (6) is zero, 
we have = uj^, which gives the critically damped 

solution: , , 

2 ;^ = cite^ ^ -h C2e^ ^ . (7) 


The columns of W are the vibrational modes of the 
object being modeled. (See figure 4.) Each mode has the 
property that a displacement or velocity over the object 
that is a scalar multiple of the mode will produce an ac- 
celeration that is also a scalar multiple of the mode. This 
property means that the modes do not interact with each 
other, which is why decoupling the system into a set of 
independent oscillators was possible. The eigenvalue for 
each mode is the ratio of the mode’s elastic stiffness to the 
mode’s mass, and it is the square of the mode’s natural 
frequency (in radians per second). In general the eigen- 
values will be positive, but for each free body in the sys- 
tem there will be six zero eigenvalues that correspond to 

^ Equivalently let W = L~^\Z where AI = LLX (Cholesky 
decomposition) and = L~^ KL~^ (symmetric eigendecom- 

position). 


the body’s six rigid-body modes. The rigid-body eigen- 
values are zero because a rigid-body displacement will 
not generate any elastic forces. 

The decoupled system of equations is not an approxi- 
mation of the original linear system, it will generate ex- 
actly the same results as the original linear system. Of 
course the linear system may have been an approxima- 
tion to some initial nonlinear one, but any problem that 
could be solved using equation (2) could also be solved 
with equation (3). Eurthermore, simulation that would 
have required numerical time integration of equation (1) 
can now be solved without integration using the analyti- 
cal solutions in equations (5) or (7). 

3.2 Discarding Modes 

Although decoupling equation (1) and then solving each 
of the resulting components analytically provides signifi- 
cant benefits, we can derive additional benefit by consid- 
ering whether or not each of these components is needed. 
In particular we can discard modes that will have no sig- 
nificant effect on the phenomena we wish to model. 

If the eigenvalue, A^, associated with a particular mode 
is large, then the force required to cause a discernible 
displacement of that mode will also be large. We can 
expect that in a given environment there will be both 
an upper bound on the magnitude of the forces encoun- 
tered and a lower limit on the amplitude of observable 
movement. Eor example, if modeling an indoor envi- 
ronment we would not expect to encounter forces in ex- 
cess of 60, 000 N (the braking force of a large truck), and 
we would not be able to observe displacements less than 
about 0.1mm. Thus if ||it?i|p/A^ < min_res/ max_frc 
for some mode then that mode’s behavior will be unob- 
servable. 

The imaginary part of uoi determines the frequency that 
a mode will vibrate at. Modes that vibrate at more than 
half the display’s frame rate will cause temporal aliasing. 

Removing modes that are too stiff and/or too high fre- 
quency to be observed will not change the appearance of 
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the resulting simulation, but removing them will greatly 
reduce the simulation’s cost. For most objects that we 
have worked with, nearly all of the modes are unobserv- 
able. A typical result is that an object with several thou- 
sand vertices will have many fewer than fifty modes that 
need to be retained. Furthermore, the number of modes 
that must be retained is nearly independent from the res- 
olution of the model. 

For later convenience let W be the matrix W with 
the columns corresponding to the discarded modes re- 
moved, and let W ^ ht the matrix W~^ with the rows 
corresponding to the discarded modes removed. Note 
that + {W)~\ W and W~ ^ are not square, and 

W~^W = I but WW~^ ^ I. 


3.3 Oscillator Coefficients and Time Steps 

The analytical solution for each mode, equation ( 5 ), de- 
scribes how that mode will behave when no external 
forces are acting on it. Using these solutions, however, re- 
quires some way of modeling responses to external forces 
and of setting initial conditions. 

Given a set of initial conditions described by the node 
positions, do, and their velocities, do, setting the oscil- 
lators to match those conditions requires finding appro- 
priate values for the coefficients ci and C2. First, the 
initial conditions are transformed to modal coordinates: 
zo = W do and zo = W do. For each mode, ci 
and C2 are given by 

( 8 ) 


{oiiXi + (^2)^0 + 2 i:o 
2^ {oLiXi + 0 ^ 2 )^ — 4A^ 


C2 


^ {c^iXi + ^2)^0 + 2 i:o 

2 2 



( 9 ) 


For the critically damped case ci and C2 are given by 

(oflAi + Q^2)^0 , . . 

Cl = + 2^0 ( 

C2 = Zo . ( 


Note that if the uf" are real then ci and C2 will also be 
real. If the ujf are complex then the ujf and the ci and 
C2 will be complex conjugate pairs. In either case equa- 
tion ( 6 ) will evaluate to a real value. 

To compute the response of a mode to an impulse de- 
livered at t = 0 , first transform the impulse to modal 
coordinates with Atg = AtW~^ f and then compute ci 
and C2 as shown above with zo set to zero and zo replaced 
by At^. Because the modes behave linearly, the response 
of the system to forces applied at an arbitrary time may 
be computed by time- shifting this impulse response and 
adding it to the existing values. 

Because ce^t-\-At)u; _ ^ 

oscillator can be stored simply as a pair of complex num- 
bers that reflect the current values of and 026^^ . 


Each time the system is advanced forward in time, these 
values get multiplied by . If At is constant then 

the step multiplier for each mode may be cached to avoid 
the cost of evaluating an exponential. Impulses applied 
to the system simply require adding the appropriate val- 
ues to each oscillator’s state. Finally, modes where 0;+ 
and uj~ are complex conjugate pairs {i.e. underdamped 
modes) can be reduced to only a single oscillator. 

3.4 Constraints 

Although we can compute the behavior of the decom- 
posed system extremely efficiently, the method is not 
particularly useful unless it accommodates manipulation 
constraints and collision response. When working with 
the original system constraints on the node positions are 
nearly trivial to implement. Collision response requires 
more sophistication but still is conceptually straightfor- 
ward. Unfortunately, applying these same constraints in 
the modal basis requires moving between the node po- 
sitions and modal coordinates which can be unintuitive. 
Matters are further complicated because if we have dis- 
carded any modes then the transformations will be non- 
invertible. 

3.4.1 Interactive Manipulation 

If we wish to include continual constraints on part of the 
system, the optimal way to do so is to remove those de- 
grees of freedom prior to performing the modal decom- 
position. Examples demonstrating this approach can be 
seen in James and Pai’s modal method for modeling tis- 
sue deformation [ 10 ], and in our deformable sheet ex- 
ample. (See accompanying animations.) Using this ap- 
proach for dynamic constraints, however, would require 
recomputing the eigendecomposition each time a con- 
straint was added or removed from the system. James 
and Pai accomplished something similar for a boundary 
element method using Sherman-Morrison- Woodbury up- 
dates but we do not know of any corresponding incremen- 
tal update scheme for an eigensystem [ 9 ]. 

Instead we apply manipulation constraints to the de- 
composed system. Let be the set of degrees of freedom 
in the original system that we wish to constrain, and let 
(j) be the places where we are willing to apply forces in 
order to enforce the constraints. For a manipulation task 
where a point on the object is being dragged we would 
typically have (p = ip but we will not require it. Let d^ 
or denote the displacement or force vectors where all 
except the elements corresponding to 0 or 7/^ have been 
removed. Similarly, let W ^ be W where all the rows 
not in have been removed and let be where 

all the columns but for those in 0 have been removed. Li- 

.. * 

nally, let d^ be the desired accelerations at the constraint 
locations. By combining id = Wz, g = f and a bit 
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of manipulation we obtain: 

dl = W^{z + wlf^). ( 12 ) 

Solving for yields: 

(ci; - W^z) , (13) 

where denotes a pseudoinverse. Velocity constraints 
only differ in that gets replaced by an impulse, e.g. 
and we have: 

f<p = Xt (4 - • ( 14 ) 

Position constraints can be enforced in a similar fashion 
so long as we adjust for how each mode will evolve over 
the interval while the force is applied: 

^ {w^swl) (d; - w^z) , (15) 

where S the diagonal matrix with components given by 

gAtw+ _ gAtw- 

+ a2)^^\ 

that compensates for the motion of each mode during the 
interval. 

3.4.2 Dynamics Simulation 

Implementing a deformable dynamics simulator for free 
bodies using modal analysis can be accomplished by 
combining the modal simulation with a standard rigid- 
body dynamics simulator. The modal system is embed- 
ded in a rigid-body reference frame, and both systems 
evolve over time. The two systems interact with each 
other though inertial effects. The modal system should 
experience centrifugal and coriolis forces as the rigid- 
body moves, and the inertial moments of the rigid-body 
will change as the modal system deforms. Unless the ob- 
ject is rotating rapidly, neither effect will be significant 
so we omit them. They could be included at an additional 
computational cost. Inertial effects due to translational 
and rotational acceleration of the rigid-body frame do not 
need to be modeled explicitly so long as the forces gen- 
erating those accelerations are also applied to the modal 
system. 

Because we are modeling deformable objects, a colli- 
sion detection method optimized for use with rigid-body 
simulations requires some modification because precom- 
puted data structures will become invalid as the object 
deforms. The method we are using employs a hierar- 
chy of axis-aligned bounding boxes, aligned to the world 
axes, to efficiently find potential collisions. The tree is 
initially constructed based on the undeformed shape of 
the object. Each leaf node in the tree corresponds to one 
of the primitives that makes up the object, and the bound- 
ing box at that node encloses the primitive. The bounding 


boxes of interior nodes encompass the union of their chil- 
dren. The tree’s topology is chosen to minimize the over- 
lap among the interior nodes. Once the object deforms 
the tree will become invalid, but recomputing the tree’s 
topology every time- step would be prohibitively expen- 
sive. Instead we use an update scheme similar to one 
described by van den Bergen [25]. After each time-step 
the bounding boxes are updated, but the tree’s topology 
does not change. If we expected arbitrary deformation, 
this could result in a very poorly structured tree, but be- 
cause the extent of deformation is limited we have found 
this approach to work quite well. 

Using these trees the collision system can efficiently 
determine contact points and a normal for each contact. 
For collisions between an object and a ground plane, the 
collision normal is simply the plane’s normal. For colli- 
sions between objects, we look at involved tetrahedra to 
determine a normal based on their overlap [16]. We have 
found that each physical contact site may produce several 
pairs of colliding primitives. To reduce the computation 
when using constraint-based collisions we cluster nearby 
collision points and treat each cluster as a single collision 
point. 

We have implemented collision response using both 
a penalty-based method and using constraints. As 
one would expect, the penalty methods require less 
work per time-step, achieving real-time performance, but 
stiff penalty coefficients can lead to instability. The 
constraint-based method requires more work per time- 
step, but it is more stable. Because the modal system will 
allow arbitrarily large time-steps in the absence of exter- 
nal infiuences we prefer the more stable constraint-based 
methods. 

To implement penalty methods, when a point on a sur- 
face violates one of the penalty constraints, a force pro- 
portional to the magnitude of the violation is applied at 
that point. Transforming the forces to modal coordinates 
and then applying the force to the modal system is done 
as described previously. The penalty force should be ap- 
plied to both the modal and the rigid-body systems. 

Constraint-based collisions require a more complex 
implementation, but we find that they produce better re- 
sults. First, when a collision occurs, the simulation is 
backed up to the point during the time- step when the ob- 
jects first came into contact. Then contact forces are cal- 
culated as the minimal outward normal force to ensure 
that the objects will not continue to penetrate. These are 
determined by solving a linear programming problem for 
the normal forces at all contact points. Baraff details an 
efficient method for solving for the required forces [1]. 

Constraint methods are often used in traditional rigid- 
body simulations only to solve for resting contact, while 
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impulses are used to calculate elastic response. Elastic 
components of the response can be handled differently in 
our modal simulation, because the elastic behavior of the 
modal system models them directly. We first enforce a 
velocity constraint that solves for an impulse to ensure 
that none of the contact velocities are negative, then sec- 
ondly it enforces an acceleration constraint that solves for 
a force to ensure that none of the contact accelerations 
are negative. The derivation of these methods requires 
equations relating the change in velocity and acceleration 
with respect to an applied impulse and acceleration, re- 
spectively. 

Let Pi be the location of a contact point on an object 
expressed in the local coordinate frame of the rigid body. 
This location will be a linear function of the modal coor- 
dinates so that: 


Pi = UWz , (17) 

where t/ is a matrix that averages the appropriate node 
locations based on the barycentric location of p in one of 
the surface triangles. The location in world coordinates 
is given by 

p^ = t + Rpi, (18) 


where t and R are the translation and rotation matrices 
for the rigid-body frame. Differentiating with respect to 
time to obtain the world velocity and acceleration of p 
yields: 

Pw = t + R[^]Pi + Rpi , (19) 

p^ = i+R[uj][u>]pi+R[a]pi+2R[u}]pi+Rpi , ( 20 ) 


where u? and a are the rigid-body’s angular velocity 
and acceleration^. The notation [a] denotes the skew- 
symmetric matrix such that [a]b = a x b = — [b]a. 

Differentiating equation (19) with respect to an applied 
impulse allows us to obtain the change in velocity gener- 
ated by a constraint force over a time interval: 


= At + R[H-Wi]pi + RUWw'^f^ 

( 21 ) 

where H is the object’s inertia matrix and r is the torque 
generated by /. Differentiating equation (20) with re- 
spect to an applied force produces a similar result for the 
change in acceleration at the contact point. These equa- 
tions are linear in /, and can be used similarly to solve 
for position, joint, and collision constraints. Position con- 
straints require that a point’s velocity and acceleration are 
zero. Joint constraints require relative velocities and ac- 
celerations are zero, merely requiring a subtraction of the 
proper terms. Collision constraints require the normal 
components of relative velocities and accelerations are 
nonnegative, and only solve for the nonnegative normal 


^In order to adhere to common convention we are reusing uj and 
a, that were previously used for the modal frequencies and Rayleigh 
damping coefficients. The intended meaning should be clear from con- 
text and the presence/absence of bold notation. 


Example 

Fig 

Verts. 

Nodes 

Tets. 

Modes 

Time 

Brain 

1 

18,847 

304 

997 

40 

68.5 sec 

Dodo 

5 

336 

113 

295 

40 

6.2 sec 

Bunny 

8 

2,633 

37,114 

15,507 

32 

24 min 

Sphere 

video 

66 

80 

282 

40 

2.9 sec 

Sheet 

video 

195 

195 

486 

20 

14.4 sec 

Bat 

video 

241 

310 

1,030 

20 

68.9 sec 


Table 1: This table list the number of vertices in the 
rendered models, the number of nodes and elements in 
the hnite element models, the number of modes retained, 
and the time required to compute the decomposition for 
some of the demonstration objects. 

force magnitude. All constraints are solved simultane- 
ously as a linear program. Solving cannot always be done 
in real-time if there are a large number of contact points, 
although system response does remain interactive. 

We model friction at the contacts using a simplified 
Coulomb friction model. The system computes a force 
opposite the tangential velocity at the contact points. The 
magnitude of the force equals the magnitude of the nor- 
mal force multiplied by a friction coefficient. If the fric- 
tion force causes the predicted tangential velocity to be 
reversed then it is limited to the force that would cause 
no slipping. If interactivity can be sacrificed, a more pre- 
cise method would be to add an additional no- slip con- 
straint to be re-solved with the other constraints. We find 
our heuristic reasonable for producing plausible friction 
effects. 

4 Results 

We have implemented a system that models deformable 
objects using a hybrid formulation that combines rigid- 
body motion with deformation computed using modal 
analysis. Objects may be interactively manipulated by 
the user with both penalty forces and displacement con- 
straints. The modal objects may collide with each other 
and with their environment. Collisions can be treated 
with either penalty forces or constraints, and objects may 
also be attached together using joint constraints. Table 1 
lists several of the models we have used to demonstrate 
our results and shows the geometric and kinematic com- 
plexity of the models along with how much precomputa- 
tion time was required to perform the modal decomposi- 
tion for each model. 

The brain model in figure 1 demonstrates pulling and 
pushing using force application. Force vectors are pro- 
jected into the modal basis, modifying the modal state, 
and then are projected out, resulting in realistic deforma- 
tion. The images in figure 6 and figure 7 show pulling and 
pushing using manipulation constraints. Typically, up to 
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Figure 5: This image sequence shows frames from an animation of a pair of objects colliding with each other. Each 
object is a hybrid simulation that incorporates a rigid and a deformable (modal) component 



Figure 6: These images shows how constraints can be 
used to deform objects. The object on the left of each im- 
age shows the object prior to deformation, and the right 
object shows the results after the red constraint points 
have been moved. 


Figure 7; These images are screen shots from an applica- 
tion running natively on a Sony PlayStation2. The yellow 
circle highlights the cursor that the user is using to poke 
and pull an elastic figure. 


around 10 points on the model can be constrained in real- 
time on a moderate speed computer (300 MHz Pentium 
II or Sony Playstation2). A limit is reached because the 
solutions to equation (13) and equation (15) require a rel- 
atively expensive computation of singular value decom- 
positions, which cannot be calculated in real-time once 
the matrices become too large. 

We have created several animations (see supplemental 
materials) demonstrating this system, each simulated in- 
teractively for moderately complex objects. The results 
appear plausible, and resemble animations that might be 
simulated using more straightforward but more compu- 
tationally expensive methods. The bottlenecks in hybrid 
modal/rigid-body simulation are collision detection and 
solving the linear program for the constraints. To reduce 
the computation used in solving the linear program, the 
extent of contact point clustering may be tweaked to sac- 
rifice accuracy for speed. Figures 5 and 8 show objects 
involved in collisions with a ground plane and each other. 


Figure 8: A sequence of images showing the Stanford 
Bunny model bouncing off a ground plane. 

As with other methods based on tetrahedral finite el- 
ements, we can embed high-resolution or non-manifold 
surfaces inside a tetrahedral volume model. The bene- 
fits of this technique are that the surface shading and tex- 
turing can be specified independently from the dynam- 
ics, and poorly constructed “polygon-soup” models may 
be used. Both the brain model in figure 1, an extremely 
complex object, and the “dodo” model in figure 5, a non- 
manifold object, are modeled in this way. The “dodo” 
model also demonstrates non-uniform material proper- 
ties: the legs and beak are made of a stiffer material than 
the rest of the body. 

5 Conclusions 

Modal analysis has been shown to be a useful tool for in- 
teractively producing realistic simulations of elastic de- 
formation. Both the analytic calculation of modal ampli- 
tudes using complex oscillators and the removal of high- 
frequency modes have a stabilizing effect on simulations, 
allowing for large time steps to be taken. 

Despite the approximation of linearity in modal anal- 
ysis, the simulation results are quite plausible for most 
objects. The exceptions are long, thin, or highly de- 
formable objects, where nonlinear behavior dominates 
the expected behavior. Despite these specific drawbacks, 
many objects can be manipulated quite efficiently and re- 
alistically using modal models. 

The already small costs of modal analysis can be re- 
duced further by leveraging graphics hardware, as shown 
by James and Pai [10] or our own implementation on the 
Sony PlayStation2. Using such hardware, CPU costs can 
be reduced to modifying mode amplitudes during evolu- 
tion of time steps, projection of forces, and application of 
manipulation constraints. 
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We recognize that there are many implementation de- 
tails that cannot fit into this paper, so we have re- 
leased the source code for our Linux implementation 
under the GNU License. It is our hope that making 
this code available will encourage others to work with 
modal simulation methods. The code may be accessed at 
www.cs.berkeley.edu/~job /Projects /ModeDef . 
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Abstract 

In this paper, we augment existing techniques for simulating flex- 
ible objects to include models for crack initiation and propagation 
in three-dimensional volumes. By analyzing the stress tensors com- 
puted over a flnite element model, the simulation determines where 
cracks should initiate and in what directions they should propagate. 
We demonstrate our results with animations of breaking bowls, 
cracking walls, and objects that fracture when they collide. By 
varying the shape of the objects, the material properties, and the 
initial conditions of the simulations, we can create strikingly dif- 
ferent effects ranging from a wall that shatters when it is hit by a 
wrecking ball to a bowl that breaks in two when it is dropped on 
edge. 

CR Categories: 1.3.5 [Computer Graphics]: Computational 

Geometry and Object Modeling — Physically based modeling; 
1.3.7 [Computer Graphics]: Three-Dimensional Graphics and 
Realism — Animation; 1.6.8 [Simulation and Modeling]: Types of 
Simulation — Animation 

Keywords: Animation techniques, physically based modeling, 

simulation, dynamics, fracture, cracking, deformation, flnite ele- 
ment method. 


1 Introduction 

With the introduction in 1998 of simulated water in Antz [5, 14] 
and clothing in GerVs Game [4, 15], passive simulation was clearly 
demonstrated to be a viable technique for commercial animation. 
The appeal of using simulation for objects without an internal 
source of energy is not surprising, as passive objects tend to have 
many degrees of freedom, making keyframing or motion capture 
difficult. Furthermore, while passive objects are often essential to 
the plot of an animation and to the appearance or mood of the piece, 
they are not characters with their concomitant requirements for con- 
trol over the subtle details of the motion. Therefore, simulations in 
which the motion is controlled only by initial conditions, physical 
equations, and material parameters are often sufficient to produce 
appealing animations of passive objects. 
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Figure 1: Slab of simulated glass that has been shattered by a heavy 
weight. 

Our approach to animating breaking objects is based on lin- 
ear elastic fracture mechanics. We model three-dimensional vol- 
umes using a flnite element method that is based on techniques 
presented in the computer graphics and mechanical engineering 
literature [3, 6, 18]. By analyzing the stresses created as a vol- 
umetric object deforms, the simulation determines where cracks 
should begin and in what directions they should propagate. The 
system accommodates arbitrary propagation directions by dynami- 
cally retesselating the mesh. Because cracks are not limited to el- 
ement boundaries, the models can form irregularly shaped shards 
and edges as they shatter. 

We demonstrate the power of this approach with the following 
examples: a glass slab that shatters when a weight is dropped onto 
it (Figure 1), an adobe wall that crumbles under the impact of a 
wrecking ball (Figure 9), a series of bowls that break when they hit 
the floor (Figure 11), and objects that break when they collide with 
each other (Figure 14). To assess the realism of this approach, we 
compare high-speed video images of a physical bowl dropping onto 
concrete and a simulated version of the same event (Figure 13). 

2 Background 

In the computer graphics literature, two previous techniques have 
been developed for modeling dynamic, deformation-induced frac- 
ture. In 1988, Terzopoulos and Fleischer [18, 19] presented a 
general technique for modeling viscoelastic and plastic deforma- 
tions. Their method used three fundamental metric tensors to de- 
flne energy functions that measured deformation over curves, sur- 
faces, and volumes. These energy functions provided the basis for 
a continuous deformation model that they simulated using a va- 
riety of discretization methods. One of their methods made use 
of a flnite differencing technique deflned by controlled continuity 
splines [17]. This formulation allowed them to demonstrate how 
certain fracture effects could be modeled by setting the elastic co- 
efficients between adjacent nodes to zero whenever the distance 
between the nodes exceeded a threshold. They demonstrated this 
technique with sheets paper and cloth that could be torn apart. 
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In 1991, Norton and his colleagues presented a technique for 
animating 3D solid objects that broke when subjected to large 
strains [12]. They simulated a teapot that shattered when dropped 
onto a table. Their technique used a spring and mass system to 
model the behavior of the object. When the distance between 
two attached mass points exceeded a threshold, the simulation sev- 
ered the spring connection between them. To avoid having flexi- 
ble strings of partially connected material hanging from the object, 
their simulation broke an entire cube of springs at once. 

Two limitations are inherent in both of these methods. First, 
when the material fails, the exact location and orientation of the 
fracture are not known. Rather the failure is deflned as the entire 
connection between two nodes, and the orientation of the fracture 
plane is left undeflned. As a result, these techniques can only re- 
alistically model effects that occur on a scale much larger than the 
inter-node spacing. 

The second limitation is that fracture surfaces are restricted to the 
boundaries in the initial mesh structure. As a result, the fracture pat- 
tern exhibits directional artifacts, similar to the “j aggies” that occur 
when rasterizing a polygonal edge. These artifacts are particularly 
noticeable when the discretization follows a regular pattern. If an ir- 
regular mesh is used, then the artifacts may be partially masked, but 
the fractures will still be forced onto a path that follows the element 
boundaries so that the object can break apart only along predeflned 
facets. 

Other relevant work in the computer graphics literature includes 
techniques for modeling static crack patterns and fractures induced 
by explosions. Hirota and colleagues described how phenomena 
such as the static crack patterns created by drying mud can be mod- 
eled using a mass and spring system attached to an immobile sub- 
strate [8]. Mazarak et al. use a voxel-based approach to model 
solid objects that break apart when they encounter a spherical blast 
wave [9]. Neff and Fiume use a recursive pattern generator to di- 
vide a planar region into polygonal shards that fly apart when acted 
on by a spherical blast wave [10]. 

Fracture has been studied more extensively in the mechanics lit- 
erature, and many techniques have been developed for simulating 
and analyzing the behavior of materials as they fail. A number of 
theories may be used to describe when and how a fracture will de- 
velop or propagate, and these theories have been employed with 
various numerical methods including flnite element and flnite dif- 
ference methods, boundary integral equations, and molecular parti- 
cle simulations. A comprehensive review of this work can be found 
in the book by Anderson [1] and the survey article by Nishioka [11]. 

Although simulation is used to model fracture both in computer 
graphics and in engineering, the requirements of the two flelds are 
very different. Engineering applications require that the simulation 
predict real-world behaviors in an accurate and reliable fashion. In 
computer animation, what matters is how the fracture looks, how 
difflcult it was to make it look that way, and how long it took. Al- 
though the technique presented in this paper was developed using 
traditional engineering tools, it is an animation technique and relies 
on a number of simpliflcations that would be unacceptable in an 
engineering context. 



Figure 2: The material coordinates deflne a 3D parameterization of 
the object. The function x{u) maps points from their location in the 
material coordinate frame to their location in the world coordinates. 
A fracture corresponds to a discontinuity in x{u). 


We derive our deformation technique by deflning a set of differ- 
ential equations that describe the aggregate behavior of the material 
in a continuous fashion, and then using a flnite element method to 
discretize these equations for computer simulation. This approach 
is fairly standard, and many different deformation models can be 
derived in this fashion. The one presented here was designed to be 
simple, fast, and suitable for fracture modeling. 

3.1 Continuous Model 

Our continuous model is based on continuum mechanics, and an ex- 
cellent introduction to this area can be found in the text by Fung [7]. 
The primary assumption in the continuum approach is that the 
scale of the effects being modeled is signiflcantly greater than the 
scale of the material’s composition. Therefore, the behavior of the 
molecules, grains, or particles that compose the material can be 
modeled as a continuous media. Although this assumption is often 
valid for modeling deformations, macroscopic fractures can be sig- 
niflcantly influenced by effects that occur at small scales where this 
assumption may not be valid. Because we are interested in graph- 
ical appearance rather than rigorous physical correctness, we will 
put this issue aside and assume that a continuum model is adequate. 

We begin the description of the continuous model by deflning 
material coordinates that parameterize the volume of space occu- 
pied by the object being modeled. Let u = [u, u, be a vector 
in that denotes a location in the material coordinate frame as 
shown in Figure 2. The deformation of the material is deflned by 
the function x{u) = [x,y,z]^ that maps locations in the material 
coordinate frame to locations in world coordinates. In areas where 
material exists, x(u) is continuous, except across a flnite number 
of surfaces within the volume that correspond to fractures in the 
material. In areas where there is no material, x(u) is undeflned. 

We make use of Green’s strain tensor, e, to measure the local 
deformation of the material [6]. It can be represented as a 3 x 3 
symmetric matrix deflned by 


€ij 


f dx dx\ 

duj J 


where 6ij is the Kronecker delta: 


( 1 ) 


3 Deformations 



( 2 ) 


Fractures arise in materials due to internal stresses created as the 
material deforms. Our goal is to model these fractures. In order 
to do so, however, we must first be able to model the deformations 
that cause them. To provide a suitable framework for modeling 
fractures, the deformation method must provide information about 
the magnitude and orientation of the internal stresses, and whether 
they are tensile or compressive. We would also like to avoid defor- 
mation methods in which directional artifacts appear in the stress 
patterns and propagate to the resulting fracture patterns. 


This strain metric only measures deformation; it is invariant with re- 
spect to rigid body transformations applied to x and vanishes when 
the material is not deformed. It has been used extensively in the 
engineering literature. Because it is a tensor, its invariants do not 
depend on the orientation of the material coordinate or world sys- 
tems. The Euclidean metric tensor used by Terzopoulos and Eleis- 
cher [18] differs only by the Sij term. 

In addition to the strain tensor, we make use of the strain rate 
tensor, i/, which measures the rate at which the strain is changing. 
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It can be derived by taking the time derivative of (1) and is defined 
by 


/ dx 

dx ^ 

\.(dx 

dx \ 

1 ^XXi 

duj j 

1 \dui 

duj J 


where an over dot indicates a derivative with respect to time sueh 
that X is the material velocity expressed in world coordinates. 

The strain and strain rate tensors provide the raw information 
that is required to compute internal elastic and damping forces, but 
they do not take into account the properties of the material. The 
stress tensor, <t, combines the basic information from the strain and 
strain rate with the material properties and determines forces inter- 
nal to the material. Like the strain and strain rate tensors, the stress 
tensor can be represented as a 3 x 3 symmetric matrix. It has two 
components: the elastic stress due to strain, and the viscous 
stress due to strain rate, The total internal stress, is the sum 

of these two components with 

O- = + o-f") . (4) 

The elastic stress and viscous stress are respectively functions of 
the strain and strain rate. In their most general linear forms, they 
are defined as 

3 3 


^ij — ^ ^ ^ ^ Cijkl €kl 

k=l 1=1 

(5) 

3 3 

k=l 1=1 

(6) 


where C is a set of the 81 elastic coefficients that relate the ele- 
ments of e to the elements and D is a. set of the 81 damping 
coefficients.^ 

Because both e and are symmetric, many of the coefficients 
in C are either redundant or constrained, and C can be reduced to 
36 independent values that relate the six independent values of e to 
the six independent values of . If we impose the additional con- 
straint that the material is isotropic, then C collapses further to only 
two independent values, /x and A, which are the Lame constants of 
the material. Equation (5) then reduces to 

3 

crlf = ^ XekkSij + 2fieij . (7) 

k=l 

The material’s rigidity is determined by the value of /x, and the 
resistance to changes in volume (dilation) is controlled by A. 

Similarly, D can be reduced to two independent values, 0 and ^ 
and ( 6 ) then reduces to 

3 

cr^j^ = ^ . ( 8 ) 

k=i 

The parameters /x and A will control how quickly the material dis- 
sipates internal kinetic energy. Since is derived from the rate 
at which e is changing, will not damp motions that are locally 
rigid, and has the desirable property of dissipating only internal vi- 
brations. 

Once we have the strain, strain rate, and stress tensors, we can 
compute the elastic potential density, 77 , and the damping potential 
density, /^, at any point in the material using 

i=l j = l 

^Actually C and D are themselves rank four tensors, and (5) and (6) are 
normally expressed in this form so that C and D will follow the standard 
rules for coordinate transforms. 



Figure 3: Given a point in the material, the traction, t, that acts 
on the surface element, dS, of a differential volume, dY, centered 
around the point with outward unit normal, n, is given by t = <t n. 

i=l j = l 

These quantities can be integrated over the volume of the material to 
obtain the total elastic and damping potentials. The elastic potential 
is the internal elastic energy of the material. The damping potential 
is related to the kinetic energy of the material after subtracting any 
rigid body motion and normalizing for the material’s density. 

The stress can also be used to compute the forces acting internal 
to the material at any loeation. Let h be an outward unit normal 
direction of a differential volume eentered about a point in the ma- 
terial. (See Figure 3.) The traction (force per unit area), t, acting 
on a face perpendicular to the normal is then given by 

t — crfi. ( 11 ) 

The examples in this paper were generated using this isotropic 
formulation. However, these techniques do not make use of the 
strain or strain rate tensors directly; rather they rely only on the 
stress. Switching to the anisotropic formulation, or even to a non- 
linear stress to strain relation, would not require any significant 
changes. 

3.2 Finite Eiement Discretization 

Before we can model a material’s behavior using this continuous 
model, it must be discretized in a way that is suitable for computer 
simulation. Two commonly used techniques are the finite difference 
and finite element methods. 

A finite difference method divides the domain of the material 
into a regular lattice and then uses numerical differencing to ap- 
proximate the spatial derivatives required to compute the strain and 
strain rate tensors. This approach is well suited for problems with 
a regular structure but becomes complieated when the structure is 
irregular. 

A finite element method partitions the domain of the material 
into distinct sub-domains, or elements as shown in Figure 4. Within 
eaeh element, the material is described loeally by a function with 
some finite number of parameters. The function is decomposed into 
a set of orthogonal shape, or basis, functions that are each associ- 
ated with one of the nodes on the boundary of the element. Adja- 
cent elements will have nodes in common, so that the mesh defines 
a piecewise function over the entire material domain. 

Our discretization employs tetrahedral finite elements with linear 
polynomial shape functions. By using a finite element method, the 
mesh can be locally aligned with the fracture surfaces, thus avoid- 
ing the previously mentioned artifacts. Just as triangles can be used 
to approximate any surface, tetrahedra can be used to approximate 
arbitrary volumes. Additionally, when tetrahedra are split along a 
fracture plane, the resulting pieces can be decomposed exactly into 
more tetrahedra. 

We chose to use linear elements because higher-order elements 
are not cost effective for modeling fracture boundaries. Although 
higher-order polynomials provide individual elements with many 
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Figure 4: Tetrahedral mesh for a simple object. In (a), only the ex- 
ternal faces of the tetrahedra are drawn; in (b) the internal structure 
is shown. 



(a) (b) 


Figure 5: A tetrahedral element is defined by its four nodes. Each 
node has (a) a location in the material coordinate system and (b) a 
position and velocity in the world coordinate system. 


degrees of freedom for deformation, they have few degrees of free- 
dom for modeling fracture because the shape of a fracture is defined 
as a boundary in material coordinates. In contrast, with linear tetra- 
hedra, each degree of freedom in the world space corresponds to a 
degree of freedom in the material coordinates. Furthermore, when- 
ever an element is created, its basis functions must be computed. 
For high-degree polynomials, this computation is relatively expen- 
sive. For systems where the mesh is constant, the cost is amortized 
over the course of the simulation. However, as fractures develop 
and parts of the object are remeshed, the computation of basis ma- 
trices can become significant. 

Each tetrahedral element is defined by four nodes. A node has 
a position in the material coordinates, m, a position in the world 
coordinates, p, and a velocity in world coordinates, v. We will refer 
to the nodes of a given element by indexing with square brackets. 
Eor example, m[ 2 ] is the position in material coordinates of the 
element’s second node. (See Eigure 5.) 

Barycentric coordinates provide a natural way to define the linear 
shape functions within an element. Let b = [&i, ^ 2 , ^ 3 , ^ 4 ]^ be 
barycentric coordinates defined in terms of the element’s material 
coordinates so that 


'u' 


■m[i] nri[2] rn[3] m[4] ' 

A. 


. 1 1 1 1 . 


These barycentric coordinates may also be used to interpolate the 
node’s world position and velocity with 


where /3 is defined by 


/3 = 


■m[i] m[2] m[3] m[4] 

.1111 


-1 


(16) 


Combining (15) with (13) and (14) yields functions that interpolate 
the world position and velocity within the element in terms of the 
material coordinates: 


x{u) = pp 1 ^^ 
x{u) = “ 

where P and V are defined as 

P = [P[llP[2]P[3]P[4]] 

V = [v[l]V[2]f[3]W[4]] • 


(17) 

(18) 

(19) 

(20) 


Note that the rows of j3 are the coefficients of the shape functions, 
and /3 needs to be computed only when an element is created or the 
material coordinates of its nodes change. Eor non-degenerate ele- 
ments, the matrix in (16) is guaranteed to be non-singular, however 
elements that are nearly co-planar will cause jS to be ill-conditioned 
and adversely affect the numerical stability of the system. 

Computing the values of e and 1 / within the element requires the 
first partials of x with respect to u: 


dx 

dui 

dx 

dui 

where 

= [(5, 


= Pl3di 

( 21 ) 

= VI3di 

( 22 ) 

;i Si2 diz 0 ]^ . 

(23) 


Because the element’s shape functions are linear, these partials are 
constant within the element. 

The element will exert elastic and damping forces on its nodes. 
The elastic force on the zth node, is defined as the negative 
partial of the elastic potential density, 77 , with respect to inte- 
grated over the volume of the element. Given /3, and the po- 
sitions in world space of the four nodes we can compute the elastic 
force by 



vol 


PjlPikC^ki 


j = l k=l 1 = 1 


(24) 


where 

vol = i [(m[ 2 ] - m[i] ) X (m[ 3 ] - m[i] )] • (m[ 4 ] - m[i] ) . (25) 

Similarly, the damping force on the zth node, is defined as 
the partial of the damping potential density, k, with respect to v 
integrated over the volume of the element. This quantity can be 
computed with 


'X 


■p[ii 

P[2] 

P[31 

P[4l’ 

A. 


. 1 

1 

1 

1 . 

x 


■V[l] 

V[2] 

V[31 

V[4]- 

A. 


. 1 

1 

1 

1 . 


(13) 

(14) 


To determine the barycentric coordinates of a point within the 
element specified by its material coordinates, we invert ( 12 ) and 
obtain 


5 = /3 


u 

1 


(15) 


= - ^ E Py] E E ■ (26) 

j=i k=i 1=1 

Summing these two forces, the total internal force that an element 
exerts on a node is 

4 33 

/w = -^Ep[iiEE IdjildikCTki 5 (27) 

j=i k=i 1=1 
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and the total internal force acting on the node is obtained by sum- 
ming the forces exerted by all elements that are attached to the node. 

As the element is compressed to less than about 30% of its ma- 
terial volume, the gradient of 77 and start to vanish causing the 
resisting forces to fall off. We have not found this to be a problem 
as even the more squishy of the materials that we have modeled 
conserve their volume to within a few percent. 

Using a lumped mass formulation, the mass contributed by an 
element to each one of its nodes is determined by integrating the 
material density, p, over the element shape function associated with 
that node. In the case of tetrahedral elements with linear shape 
functions, this mass contribution is simply p vol/4. 

The derivations above are sufficient for a simulation that uses an 
explicit integration scheme. Additional work, including computing 
the Jacobian of the internal forces, is necessary for implicit integra- 
tion scheme. (See for example [2] and [3].) 

3.3 Collisions 

In addition to the forces internal to the material, the system com- 
putes collision forces. The collision forces are computed using a 
penalty method that is applied when two elements intersect or if an 
element violates another constraint such as the ground. Although 
penalty methods are often criticized for creating stiff equations, we 
have found that for the materials we are modeling the internal forces 
are at least as stiff as the penalty forces. Penalty forces have the 
advantage of being very fast to compute. We have experimented 
with two different penalty criteria: node penetration and overlap 
volume. The examples presented in this paper were computed with 
the node penetration criteria; additional examples on the conference 
proceedings CD-ROM were computed with the overlap volume cri- 
teria. 

The node penetration criteria sets the penalty force to be pro- 
portional to the distance that a node has penetrated into another 
element. The penalty force acts in the direction normal to the pene- 
trated surface. The reaction force is distributed over the penetrated 
element’s nodes so that the net force and moment are the negation of 
the penalty force and moment acting at the penetrating node. This 
test will not catch all collisions, and undetected intersecting tetra- 
hedra may become locked together. It is however, fast to compute, 
easy to implement, and adequate for situations that do not involve 
complex collision interactions. 

The overlap volume criteria is more robust than the node pene- 
tration method, but it is also slower to compute and more complex 
to implement. The intersection of two tetrahedral elements is com- 
puted by clipping the faces of each tetrahedron against the other. 
The resulting polyhedron is then used to compute the volume and 
center of mass of the intersecting region. The area weighted nor- 
mals of the faces of the polyhedron that are contributed by one of 
the tetrahedra are summed to compute the direction that the penalty 
force acts in. A similar computation can be performed for the other 
tetrahedra, or equivalently the direction can be negated. Provided 
that neither tetrahedra is completely contained within the other, this 
criteria is more robust than the node penetration criteria. Addition- 
ally, the forces computed with this method do not depend on the 
object tessellation. 

Computing the intersections within the mesh can be very expen- 
sive, and we use a bounding hierarchy scheme with cached traver- 
sals to help reduce this cost. 


4 Fracture Modeling 

Our fracture model is based on the theory of linear elastic fracture 
mechanics [1]. The primary distinction between this and other the- 
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Figure 6 : Three loading modes that can be experienced by a crack. 
Mode I: Opening, Mode II: In-Plane Shear, and Mode III: Out-of- 
Plane Shear. Adapted from Anderson [1]. 


ories of fracture is that the region of plasticity near the crack tip^ is 
neglected. Because we are not modeling the energy dissipated by 
this plastic region, modeled materials will be brittle. This statement 
does not mean that they are weak; rather the term “brittle” refers to 
the fact that once the material has begun to fail, the fractures will 
have a strong tendency to propagate across the material as they are 
driven by the internally stored elastic energy. 

There are three loading modes by which forces can be applied to 
a crack causing it to open further. (See Figure 6 .) In most circum- 
stances, some combination of these modes will be active, producing 
a mixed mode load at the crack tip. For all three cases, as well as 
mixed mode situations, the behavior of the crack can be resolved 
by analyzing the forces acting at the crack tip: tensile forces that 
are opposed by other tensile forces will cause the crack to continue 
in a direction that is perpendicular to the direction of largest tensile 
load, and conversely, compressive loads will tend to arrest a crack 
to which they are perpendicular. 

The finite element model describes the surface of a fracture with 
elements that are adjacent in material coordinates but that do not 
share nodes across the internal surface. The curve that represents 
the crack tip is then implicitly defined in a piecewise linear fashion 
by the nodes that border the fracture surface, and further extension 
of the crack may be determined by analyzing the internal forces 
acting on these nodes. 

We will also use the element nodes to determine where a crack 
should be initiated. While this strategy could potentially introduce 
unpleasant artifacts, we note that because the surface of an object is 
defined by a polygonal boundary (the outer faces of the tetrahedra) 
there will always be a node located at any concavities. Because con- 
cavities are precisely the locations where cracks commonly begin, 
we believe that this decision is acceptable. 

Our fracture algorithm is as follows: after every time step, the 
system resolves the internal forces acting on all nodes into their ten- 
sile and compressive components, discarding any unbalanced por- 
tions. At each node, the resulting forces are then used to form a ten- 
sor that describes how the internal forces are acting to separate that 
node. If the action is sufficiently large, the node is split into two dis- 
tinct nodes and a fracture plane is computed. All elements attached 
to the node are divided along the plane with the resulting tetrahe- 
dra assigned to one or the other incarnations of the split node, thus 
creating a discontinuity in the material. Any cached values, such as 
the node mass or the element shape functions, are recomputed for 
the affected elements and nodes. With this technique, the location 
of a fracture or crack tip need not be explicitly recorded unless this 
information is useful for some other purpose, such as rendering. 


^The term “crack tip” implies that the fracture will have a single point 
at its tip. In general, the front of the crack will not be a single point; rather 
it will be a curve that defines the boundary of the surface discontinuity in 
material coordinates. (See Figure 4.) Nevertheless, we will refer to this 
front as the crack tip. 
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4.1 Force Decomposition 

The forces acting on a node are decomposed by first separating 
the element stress tensors into tensile and compressive components. 
For a given element in the mesh, let v* (<t) , with i G { 1 , 2, 3}, be the 
zth eigenvalue of <t, and let n*(<T) be the corresponding unit length 
eigenvector. Positive eigenvalues correspond to tensile stresses and 
negative ones to compressive stresses. Since cr is real and symmet- 
ric, it will have three real, not necessarily unique, eigenvalues. In 
the case where an eigenvalue has multiplicity greater than one, the 
eigenvectors are selected arbitrarily to orthogonally span the appro- 
priate subspace [13]. 

Given a vector a in 3^^^, we can construct a 3 x 3 symmetric ma- 
trix, m(a) that has \a\ as an eigenvalue with a as the corresponding 
eigenvector, and with the other two eigenvalues equal to zero. This 
matrix is defined by 

= ^ ( 28 ) 
^ ^ ] 0 : a = 0 . 

The tensile component, <t^, and compressive component, <t“, 
of the stress within the element can now be computed by 

3 

max(0 , v^(<t)) m(n*(<T)) (29) 

i=l 

3 

<T~ = min(Q, v^(<t)) m(ti*(<T)) . (30) 

i=l 

Using this decomposition, the force that an element exerts on a 
node can be separated into a tensile component, and a com- 
pressive component, /JT^. This separation is done by reevaluating 

the internal forces exerted on the nodes using (27) with or <t“ 
substituted for <t. Thus the tensile component is 

ft] = Y.P]i] • ( 31 ) 

j=l k=l 1=1 

The compressive component can be computed similarly, but be- 
cause <T = <T^+<T“,it can be computed more efficiently using 
f [i] = /^ + / [i]- 

Each node will now have a set of tensile and a set of compressive 
forces that are exerted by the elements attached to it. For a given 
node, we denote these sets as {/^} and {f~} respectively. The 
unbalanced tensile load, is simply the sum over {/^}, and the 
unbalanced compressive load, is the sum over {f~}. 

4.2 The Separation Tensor 

We describe the forces acting at the nodes using a stress variant that 
we call the separation tensor, The separation tensor is formed 
from the balanced tensile and compressive forces acting at each 
node and is computed by 

(-m(/+)+y]m(/) + m(/-)-y^m(/) j . (32) 

\ fe{f+} fe{f-} J 

It does not respond to unbalanced actions that would produce a rigid 
translation, and is invariant with respect to transformations of both 
the material and world coordinate systems. 

The separation tensor is used directly to determine whether a 
fracture should occur at a node. Let be the largest positive eigen- 
value of If is greater than the material toughness, r, then the 



Figure 7: Diagram showing how an element is split by the fracture 
plane, (a) The initial tetrahedral element, (b) The splitting node 
and fracture plane are shown in blue, (c) The element is split along 
the fracture plane into two polyhedra that are then decomposed into 
tetrahedra. Note that the two nodes created from the splitting node 
are co-located, the geometric displacement shown in (c) only illus- 
trates the location of the fracture discontinuity. 



Figure 8: Elements that are adjacent to an element that has been 
split by a fracture plane must also be split to maintain mesh consis- 
tency. (a) Neighboring tetrahedra prior to split, (b) Face neighbor 
after split, (c) Edge neighbor after split. 


material will fail at the node. The orientation in world coordinates 
of the fracture plane is perpendicular to n+, the eigenvector of ^ 
that corresponds to . In the case where multiple eigenvalues are 
greater than r, multiple fracture planes may be generated by first 
generating the plane for the largest value, remeshing (see below), 
and then recomputing the new value for ^ and proceeding as above. 


4.3 Local Remeshing 

Once the simulation has determined the location and orientation of 
a new fracture plane, the mesh must be modified to reflect the new 
discontinuity. It is important that the orientation of the fracture be 
preserved, as approximating it with the existing element boundaries 
would create undesirable artifacts. To avoid this potential difficulty, 
the algorithm remeshes the local area surrounding the new fracture 
by splitting elements that intersect the fracture plane and modifying 
neighboring elements to ensure that the mesh stays self-consistent. 

First, the node where the fracture originates is replicated so that 
there are now two nodes, and n~ with the same material posi- 
tion, world position, and velocity. The masses will be recalculated 
later. The discontinuity passes “between” the two co-located nodes. 
The positive side of the fracture plane is associated with and the 
negative side with n~ . 

Next, all elements that were attached to the original node are ex- 
amined, comparing the world location of their nodes to the fracture 
plane. If an element is not intersected by the fracture plane, then 
it is reassigned to either orn“ depending on which side of the 
plane it lies. 

If the element is intersected by the fracture plane, it is split along 
the plane. (See Figure 7.) A new node is created along each edge 
that intersects the plane. Because all elements must be tetrahedra, in 
general each intersected element will be split into three tetrahedra. 
One of the tetrahedra will be assigned to one side of the plane and 
the other two to the other side. Because the two tetrahedra that 
are on the same side of the plane both share either or n~ , the 
discontinuity does not pass between them. 

In addition to the elements that were attached to the original 
node, it may be necessary to split other elements so that the mesh 


142 


Computer Graphics Proceedings, Annual Conference Series, 1999 



Figure 9: Two adobe walls that are struck by wrecking balls. Both walls are attached to the ground. The ball in the second row has 50 x 
the mass of the first. Images are spaced apart 133.3 ms in the first row and 66.6 ms in the second. The rightmost images show the final 
configurations. 



Figure 10: Mesh for adobe wall, (a) The facing surface of the initial 
mesh used to generate the wall shown in Figure 9. (b) The mesh af- 
ter being struck by the wrecking ball, reassembled, (c) Same as (b), 
with the cracks emphasized, (d) Internal fractures shown as wire- 
frame. 


stays consistent. In particular, an element must be split if the face or 
edge between it and another element that was attached to the orig- 
inal node has been split. (See Figure 8.) To prevent the remeshing 
from cascading across the entire mesh, these splits are done so that 
the new tetrahedra use only the original nodes and the nodes cre- 
ated by the intersection splits. Because no new nodes are created, 
the effect of the local remeshing is limited to the elements that are 
attached to the node where the fracture originated and their imme- 
diate neighbors. Because the tetrahedra formed by the secondary 
splits do not attach to either or n“, the discontinuity does not 
pass between them. 

Finally, after the local remeshing has been completed, any 
cached values that have become invalid must be recomputed. In 
our implementation, these values include the element basis matrix, 
/3, and the node masses. 

Two additional subtleties must also be considered. The first 
subtlety occurs when an intersection split involves an edge that 
is formed only by tetrahedra attached to the node where the crack 
originated. When this happens, the fracture has reached a boundary 
in the material, and the discontinuity should pass through the edge. 
Remeshing occurs as above, except that two nodes are created on 
the edge and one is assigned to each side of the discontinuity. 

Second, the fracture plane may pass arbitrarily close to an exist- 
ing node producing arbitrarily ill-conditioned tetrahedra. To avoid 
this, we employ two thresholds, one the distance between the frac- 


ture plane and an existing node, and the other on the angle between 
the fracture plane and a line from the node where the split origi- 
nated to the existing node. If either of these thresholds are not met, 
then the intersection split is snapped to the existing node. In our 
results, we have used thresholds of 5 mm and 0.1 radians. 


5 Results and Discussion 

To demonstrate some of the effects that can be generated with this 
fracture technique, we have animated a number of scenes that in- 
volve objects breaking. Figure 1 shows a plate of glass that has had 
a heavy weight dropped on it. The area in the immediate vicinity 
of the impact has been crushed into many small fragments. Further 
away from the weight, a pattern of radial cracks has developed. 

Figure 9 shows two walls being struck by wrecking balls. In 
the first sequence, the wall develops a network of cracks as it ab- 
sorbs most of the ball’s energy during the initial impact. In the sec- 
ond sequence, the ball’s mass is 50 x greater, and the wall shatters 
when it is struck. The mesh used to generate the wall sequences is 
shown in Figure 10. The initial mesh contains only 338 nodes and 
1109 elements. By the end of the sequence, the mesh has grown 
to 6892 nodes and 8275 elements. These additional nodes and el- 
ements are created where fractures occur; a uniform mesh would 
require many times this number of nodes and elements to achieve a 
similar result. 

Figure 1 1 shows the final frames from four animations of bowls 
that were dropped onto a hard surface. Other than the toughness, 
T, of the material, the four simulations are identical. The first bowl 
develops only a few cracks; the weakest breaks into many pieces. 

Because this system works with solid tetrahedral volumes rather 
than with the polygonal boundary representations created by most 
modeling packages, boundary models must be converted before 
they can be used. A number of systems are available for creating 
tetrahedral meshes from polygonal boundaries. The models that 
we used in these examples were generated either from a CSG de- 
scription or a polygonal boundary representation using NETGEN, 
a publicly available mesh generation package [16]. 

Although our approach avoids the “jaggy” artifacts in the frac- 
ture patterns caused by the underlying mesh, there remain ways in 
which the results of a simulation are infiuenced by the mesh struc- 
ture. The most obvious is that the deformation of the material is 
limited by the degrees of freedom in the mesh, which in turn limits 
how the material can fracture. This limitation will occur with any 
discrete system. The technique also limits where a fracture may ini- 
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Figure 11: Bowls with successively lower toughness values, r. Each of the bowls were dropped from the same height. Other than r, the 
bowls have same material properties. 



Figure 12: Back-cracking during fracture advance. The dashed line 
is the axis of the existing crack. Cracks advance by splitting ele- 
ments along a fracture plane, shown as a solid line, computed from 
the separation tensor, (a) If the crack does not turn sharply, then 
only elements in front of the tip will be split, (b) If the crack turns 
at too sharp an angle, then the backwards direction may not fall 
inside of the existing failure and a spurious bifurcation will occur. 


date by examining only the existing nodes. This assumption means 
that very coarse mesh sizes might behave in an unintuitive fash- 
ion. However, nodes correspond to the locations where a fracture 
is most likely to begin; therefore, with a reasonable grid size, this 
limitation is not a serious handicap. 

A more serious limitation is related to the speed at which a crack 
propagates. Currently, the distance that a fracture may travel during 
a time step is determined by the size of the existing mesh elements. 
The crack may either split an element or not; it cannot travel only 
a fraction of the distance across an element. If a crack were being 
opened slowly by an applied load on a model with a coarse resolu- 
tion mesh, this limitation would lead to a “button popping” effect 
where the crack would travel across one element, pause until the 
stress built up again, and then move across the next element. A 
second type of artifact may occur if the crack’s speed should be 
significantly greater than the element width divided by the simula- 
tion time step. In this case, a high stress area will race ahead of 
the crack tip, causing spontaneous failures to occur in the material. 
Although we have not observed these phenomena in our examples, 
developing an algorithm that allows a fracture to propagate arbitrary 
distances is an area for future work. 

Another limitation stems from the fact that while the fracture 
plane’s orientation is well defined, the crack tip’s forward direction 
is not. As shown in Figure 12, if the cracks turns at an angle greater 
than half the angle at the crack tip, then a secondary fracture will 
develop in the opposite direction to the crack’s advance. While this 
effect is likely present in some of our examples, it does not appear 
to have a significant impact on the quality of the results. If the arti- 
facts were to be a problem, they could be suppressed by explicitly 
tracking the fracture propagation directions within the mesh. 

The simulation parameters used to generate the examples in this 
paper are listed in Table 1 along with the computation time required 
to generate one second of animation. While the material density 
values, p, are appropriate for glass, stone, or ceramic, we used val- 
ues for the Lame constants, A and /i, that are significantly less than 
those of real materials. Larger values would make the simulated 
materials appear stiff er, but would also require smaller time steps. 


The values that we have selected represent a compromise between 
realistic behavior and reasonable computation time. 

Our current implementation can switch between either a forward 
Euler integration scheme or a second order Taylor integrator. Both 
of these techniques are explicit integration schemes, and subject to 
stability limits that require very small time steps for stiff materi- 
als. Although semi-implicit integration methods have error bounds 
similar to those of explicit methods, the semi-implicit integrators 
tend to drive errors towards zero rather than infinity so that they are 
stable at much larger time steps. Other researchers have shown that 
by taking advantage of this property, a semi-implicit integrator can 
be used to realize speed ups of two or three orders of magnitude 
when modeling object deformation [2]. Unfortunately, it may be 
difficult to realize these same improvements when fracture prop- 
agation is part of the simulation. As discussed above, the crack 
speed is limited in inverse proportion to the time step size, and the 
large time steps that might be afforded by a semi-implicit integrator 
could cause spontaneous material failure to proceed crack advance. 
We are currently investigating how our methods may be modified 
to be compatible with large time step integration schemes. 

Many materials and objects in the real world are not homoge- 
neous, and it would be interesting to develop graphical models for 
animating them as they fail. For example, a brick wall is made 
up of mortar and bricks arranged in a regular fashion, and if simu- 
lated in a situation like our wall example, a distinct pattern would 
be created. Similarly, the connection between a handmade cup and 
its handle is often weak because of the way in which the handle is 
attached. 

One way to assess the realism of an animation technique is by 
comparing it with the real world. Figure 13 shows high-speed video 
footage of a physical bowl as it falls onto its edge compared to our 
imitation of the real-world scene. Although the two sets of fracture 
patterns are clearly different, the simulated bowl has some qualita- 
tive similarities to the real one. Both initially fail along the leading 
edge where they strike the ground, and subsequently develop verti- 
cal cracks before breaking into several large pieces. 
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Example 

Eigure 

A (AT/m^) 


Material Parameters 
(j){Ns/rrF) ipiNs/m^) 

p (kg/m^) 

r {N/rrF) 

Minu 
Time p( 
Minimum 

tes of Computa 
2 r Simulation Si 

Maximum 

tion 

econd 

Average 

Glass 

1 

1.04 X 10® 

1.04 X 10® 

0 

6760 

2588 

10140 

75 

667 

273 

Wall #1 

9.a 

6.03 X 10® 

1.21 X 10® 

3015 

6030 

2309 

6030 

75 

562 

399 

Wall #2 

9.b 

0 

1.81 X 10® 

0 

18090 

2309 

6030 

75 

2317 

1098 

Bowl #1 

ll.a 

2.65 X 10® 

3.97 X 10® 

264 

397 

1013 

52.9 

90 

120 

109 

Bowl #2 

ll.b 

2.65 X 10® 

3.97 X 10® 

264 

397 

1013 

39.6 

82 

135 

115 

Bowl #3 

ll.c 

2.65 X 10® 

3.97 X 10® 

264 

397 

1013 

33.1 

90 

150 

127 

Bowl #4 

ll.d 

2.65 X 10® 

3.97 X 10® 

264 

397 

1013 

13.2 

82 

187 

156 

Comp. Bowl 

13 

0 

5.29 X lO"^ 

0 

198 

1013 

106 

247 

390 

347 

The End 

14 

0 

9.21 X 10® 

0 

9.2 

705 

73.6 

622 

6667 

4665 


Table 1: Material parameters and simulation times for examples. The times listed reflect the total number of minutes required to compute one 
second of simulated data, including graphics and file I/O. Times were measured on an SGI 02 with a 195 MHz MIPS RlOK processor. 



Figure 13: Comparison of a real-world event and simulation. The top row shows high-speed video images of a physical ceramic bowl dropped 
from approximately one meter onto a hard surface. The bottom row is the output from a simulation where we attempted to match the initial 
conditions of the physical bowl. Video images are 8 ms apart. Simulation images are 13 ms apart. 
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Abstract 

In this paper, we describe a method for realistically animating duc- 
tile fracture in common solid materials such as plastics and metals. 
The effects that characterize ductile fracture occur due to interac- 
tion between plastic yielding and the fracture process. By modeling 
this interaction, our ductile fracture method can generate realistic 
motion for a much wider range of materials than could be real- 
ized with a purely brittle model. This method directly extends our 
prior work on brittle fracture [O’Brien and Hodgins, SIGGRAPH 
99]. We show that adapting that method to ductile as well as brit- 
tle materials requires only a simple to implement modification that 
is computationally inexpensive. This paper describes this modifi- 
cation and presents results demonstrating some of the effects that 
may be realized with it. 

CR Categories: 1.3.5 [Computer Graphics]: Computational 

Geometry and Object Modeling — Physically based modeling; 
1.3.7 [Computer Graphics]: Three-Dimensional Graphics and 
Realism — Animation; 1.6.8 [Simulation and Modeling]: Types of 
Simulation — Animation 

Keywords: Animation techniques, physically based modeling, 

simulation, dynamics, fracture, cracking, deformation, finite ele- 
ment method, ductile fracture, plasticity. 

1 Introduction 

As techniques for generating photorealistic computer rendered im- 
ages have improved, the use of physically based animation to gen- 
erate special effects in film, television, and games has become in- 
creasingly common. Physically based animation techniques have 
proven to be particularly useful for violent or destructive effects 
that would be impractical or expensive to achieve using other meth- 
ods. For example, when creating effects for the film Pearl Har- 
bor, Industrial Light and Magic made extensive use of simulation 
methods for modeling the destruction of ships, planes, and other 
structures [Duncan, 2001]. 

Animating objects as they break, crack, tear, or in general frac- 
ture appears to be an obvious place where physically based mod- 
eling should be useful, particularly if the object is expensive, ir- 
replaceable, or if breaking it would be hazardous. However even 
the most general of current techniques for animating fracture are 
limited to modeling only brittle materials. 

The term brittle does not mean that a material is fragile. It means 
that the material experiences only elastic deformation before frac- 
ture. Few real materials are truly brittle. In contrast, ductile ma- 
terials behave elastically up to a point and then experience some 
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Figure 1 : Four hollow balls that have been dropped onto a hard 
surface. The ball in (a) flattens out and visibly demonstrates plastic 
yielding. The other three do not show an appreciable amount of 
plastic deformation, but the manner in which they split and tear, as 
opposed to shattering, arises from of the interaction between plastic 
yielding and the fracture process. 


amount of plastic deformation before fracture. When brittle materi- 
als fracture, they shatter. However, ductile materials demonstrate a 
much wider range of fracture behaviors. (See figures 1 and 3.) This 
wider range of behaviors arises due to the interaction of plastic en- 
ergy absorption with the fracture process. 

This paper describes a method suitable for modeling ductile frac- 
ture in common solid materials such as plastics or metals. This 
method directly extends our prior technique presented in [O’Brien 
and Hodgins, 1999] for modeling brittle fracture. Adapting that 
technique to ductile as well as brittle materials requires only a sim- 
ple to implement and computationally inexpensive modification. 
This extension dramatically expands the range of materials that may 
be modeled. For the sake of brevity, this paper describes only this 
modification and presents results demonstrating some of the effects 
that may be realized with it. 

2 Related Work 

The primary contribution of this paper is extending [O’Brien and 
Hodgins, 1999] to include ductile fracture by adding a plastic- 
ity model to the underlying finite-element method. The plastic- 
ity model we describe is not novel. It consists of the von Mises 
yield criterion, simple kinematic work hardening, and a finite yield 
limit [Fung, 1965]. This plasticity model is similar to the one used 
in [Terzopoulos and Fleischer, 1988a] and [Terzopoulos and Fleis- 
cher, 1988b]. The primary differences between their model and the 
one presented here are that this model realistically preserves vol- 
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ume and it includes a second elastic regime once a limit on the 
amount of plastic flow has been exceeded. 

Although ductile fracture has not been widely addressed in the 
graphics literature, several other graphics researchers have investi- 
gated brittle fracture. In [Terzopoulos and Fleischer, 1988a] and 
[Terzopoulos and Fleischer, 1988b] a finite differencing scheme 
was used to model tearing sheets of cloth-like material. Work 
by [Norton et al., 1991] used a mass/spring system to model a 
breaking teapot. Fracture in the eontext of explosions was explored 
by [Mazarak et al., 1999], [Neff and Fiume, 1999], and [Yngve 
et al., 2000]. Most recently [Smith et al., 2001] used eonstraint- 
based methods for modeling brittle fracture. 

Outside the graphics literature, both brittle and duetile fracture 
have been investigated extensively. A comprehensive review of this 
work can be found in [Anderson, 1995]. 

3 Ductile vs. Brittle Fracture 

The common usage of the terms elastic and brittle differs substan- 
tially from their teehnical meanings. For example, elastic is often 
used incorrectly as a synonym fox flexible, and the term brittle as a 
synonym for /rag//^ [Merriam- Webster, 1998]. The technically cor- 
rect definition of an elastic material refers to a material that returns 
to its original configuration when deforming forces have been re- 
moved. The ratio between the magnitude of a force and the amount 
of deformation it induees, that is how easily the material deforms, is 
the compliance of the material and it is irrelevant to whether or not 
the material is elastic. Although no real material is perfectly elas- 
tic, both natural rubber and common glass are examples of nearly 
elastic materials. Rubber’s elastic behavior is obvious while glass 
appears to be rigid. A brittle material is simply one that behaves 
elastically up until the point where it fractures. 

In contrast to an elastic material, a plastic material will not re- 
turn to its original configuration once deforming forees have been 
removed. When a material, sueh as lead, bends and then holds its 
new shape, it demonstrates plastic behavior. As previously stated, 
real materials do not behave perfectly elastically. Real materials 
can be deformed only to a limited extent before they will no longer 
return to their original configuration. This limit is known as the ma- 
terial’s elastic limit or yield point. When the elastic limit has been 
exceeded, the material enters a plastic regime and begins to experi- 
ence plastic fiow. Eventually, at the failure threshold, it fractures. 

The terms brittle and ductile relate to the relative values of the 
elastic limit and failure threshold. If the failure threshold nearly 
coincides with the elastic limit, then the material will experience 
only negligible plastic deformation before fracture. The term brit- 
tle refers to such a material. In contrast, for a ductile material the 
failure threshold is signifieantly larger than the elastic limit so that 
as the material deforms it experiences an elastic regime, followed a 
plastic regime, and then finally fracture. 

The significance of the distinction between duetile and brittle 
materials arises because elastic deformation stores energy whereas 
plastic deformation dissipates it. When a brittle material is de- 
formed to its failure threshold, the majority of the energy used to 
deform it has been stored as elastic potential. When fracture oc- 
curs, the energy is released and it tends to drive the fracture further 
into the material. Thus, even though a large or small force may 
be required to start a crack in a brittle material (depending on its 
toughness), once the crack is started only a small amount of energy 
is required to push it further. In contrast, a ductile material requires 
significantly more work to propagate a crack because energy is be- 
ing absorbed by plastic deformation. As a result, brittle materials 
tend to shatter, whereas duetile ones tend to tear. 

In general the underlying causes of plasticity are fairly compli- 
cated and they give rise to a number of phenomena. For example, 
the energy absorbed by plastie deformation does not simply vanish 
and it may result in effeets such as fatigue weakening. However for 
the purposes of animating failure events that oecur over relatively 
short periods of time, the most significant effect of plasticity is how 
it directly effects fracture propagation, and the methods discussed 
here focus on modeling those effects efficiently. Additional infor- 
mation about mathematical models of deformation and plasticity 
can be found in [Fung, 1965; Fung, 1969] and [Han and Reddy, 


1999]. Additional information concerning both brittle and ductile 
fracture may be found in [Anderson, 1995]. 

4 Modeling Ductility 

The dynamic fracture propagation teehnique described in [O’Brien 
and Hodgins, 1999] models the fracture proeess using a simple 
tetrahedral finite-element method, rules for fracture initiation and 
propagation, and procedures for automatic remeshing as a erack 
advanees. The quality of the results produced with that method 
is sufficient for graphical applications, and the only limitation that 
makes it unsuitable for modeling fracture in ductile materials is that 
the continuum model does not account for plastic deformation. 

Extending that model to account for plasticity may be aceom- 
plished by simply redefining the strain metric used to compute ele- 
ment stresses. This change has only a local impact on the fracture 
algorithm, and so we will not repeat the details of the method which 
appear in [O’Brien and Hodgins, 1999]. Instead we describe only 
the modifications that should be made to the algorithm: 

• The elastic strain, e® defined in section 4.1 of this paper, takes 
the place of the total strain, e, when computing the elastic 
stress. 

• A routine for updating the plastic strain, described in sec- 
tion 4.2 of this paper, must be called during every integration 
step. 

Even though this extension requires only incremental modifications 
to the previous method, it signifieantly extends the range of materi- 
als that may be realistically modeled. Furthermore, as our examples 
will demonstrate, small amounts of plastic yielding can dramati- 
cally effect the overall appearance of fracture patterns in a material, 
even though the plastic deformation itself cannot be observed di- 
rectly. We feel that the significant relationship between plasticity 
and the appearance of fracture in most materials makes modeling 
plasticity a required component of any general system for animat- 
ing fracture. 

4.1 Decomposing Strain 

The first step towards modeling plastie deformation requires sepa- 
rating the strain into two components: 

e = e^ + e^ (1) 

where e is the total strain, is the strain due to plastic defor- 
mation, and e® is the strain due to elastic deformation. The total 
strain is a purely geometric measure, it indicates how much the 
local shape of an object has changed from some initial referenee 
configuration and it may be computed from the material’s current 
configuration. (See [O’Brien and Hodgins, 1999] for computation 
of Green’s strain.) The plastic strain refiects how the material’s rest 
shape has been permanently distorted and it is part of the material’s 
state. Initially, the plastic strain is zero^ and it will evolve accord- 
ing to an update rule as the simulation progresses. Because the total 
and plastic strains are known at any given time, equation (1) may 
be used to compute the elastic strain. 

4.2 Plastic Update 

The algorithm for modeling the evolution of the plastic strain con- 
sists of a yield condition that must be met before plastic deforma- 
tion occurs and a rule for computing plastic fiow once the yield 
criterion has been met. We employ von Mises’s yield criterion for 
the eondition under which plastic flow will begin [Fung, 1965]. Our 
method for updating the plastic strain assumes that the rate of plas- 
tic fiow in the material is close enough to its rate of deformation so 
that plastic fiow can be updated instantaneously. This assumption 
precludes modeling phenomena such as creep and relaxation, but 

^ A non-zero initial value for the plastic strain could be used to model 
an object that has already experienced plastic deformation. 


292 



Computer Graphics Proceedings, Annual Conference Series, 2002 




Plastic Limit, ^2 
Elastic limit, 

Zero strain point 

P 

Plastic strain, 6 

Current total strain deviation, 6 


+ eP 


Figure 2 : These diagrams illustrate the behavior of the plasticity 
model, (a) Elastic deformation, (b) and (c) Plastic deformation, 
(d) Limit of plastic yield. (See explanation in the text.) 


under most circumstances these phenomena do not significantly ef- 
fect fracture behavior. We also ignore the weakening of a material 
due to repeated plastic deformation known as fatigue. While fa- 
tigue often plays a significant role in the failure of mechanisms and 
structures, a previously fatigued object may be modeled by locally 
adjusting its toughness and plastic limits. 

The von Mises yield criterion is based on the deviation of the 
elastic strain given by 


e = e — 




( 2 ) 


where Tr (•) is the trace of a matrix and I is the identity matrix. 
By averaging out the sum of the diagonal terms, the elastic strain 
deviation refiects only the portion of the elastic strain that is due to 
shape distortion and it excludes dilation. Excluding dilation makes 
the plastic deformation insensitive to hydrostatic pressure and will 
prevent the material from changing its volume which would gener- 
ate unnatural behavior. 

The yield criterion compares the magnitude of the elastic strain 
deviation (Erobenius norm) to a material constant, 71 : 


71 < 1711 . 


( 3 ) 


Together equations ( 2 ) and ( 3 ) define the von Mises yield crite- 
rion [Eung, 1965 ]. If this condition is met then plastic deformation 
will occur. We compute the base change in plastic deformation ac- 
cording to: 




He'll -71 / 
He'll 


( 4 ) 


A limit on the total amount of plastic deformation that can be with- 
stood by the material, 72, is enforced by updating the plastic strain 
at every time- step according to: 


:= (e^ + Ae^) min ( 1 


72 


\\eP-^AeP\ 


( 5 ) 


The behavior of this plasticity model is illustrated by figure 2 . 
The image plane represents a two-dimensional projection of the 
five-dimensional space of strain deviations.^ The plastic strain be- 
haves as if it were being dragged by the total strain using a rope of 
length 71. The difference between the plastic strain’s and the total 

^ For three-dimensional objects, strain is a 3 x 3 symmetric tensor with 
nine components. Because of symmetry, only six of these components are 
independent. Equation (2) removes one degree of freedom, leaving five. 



Figure 3 : Images showing the results of simulating a set of eight 
thin walls with different material parameters as they are each struck 
by a heavy projectile. A purely brittle material is shown in the top- 
left. The others images demonstrate how varying the plasticity of 
the material can produce a range of effects. 


strain’s locations represents the current elastic strain. A barrier at 
radius 72 restricts the motion of the plastic strain, but not the total 
strain. An elastic force (stress) attracts the total strain to the plas- 
tic strain, but not the plastic strain to the total strain. As shown in 
figure 2.C, the plastic deformation will depend on the history of the 
total strain’s movement. 

5 Results and Discussion 

Figure 3 shows a set of thin walls that have been struck by a heavy 
weight. The walls are clamped at the bottom, and they experi- 
ence collision forces due to contacts with the ground plane, the 
weight, and self-collisions. The top-left image in figure 3 with 
(71 = 72 = 0 ) shows the behavior of a purely brittle material. The 
other images in figure 3 show some examples that demonstrate the 
effects of different plastic parameter values. In the left column 71 
has been varied while 72 was held fixed. The right column demon- 
strates the result of varying 72 while 71 was held fixed. Some of the 
images, such as the bottom-right with (71 = 0 . 001 , 72 = 0 . 486 ), 
demonstrate obvious amounts of plastic yielding. However, plastic- 
ity also plays a significant role in the images where plastic yielding 
is not obviously visible. For example, (71 = 0.001,72 = 0 . 162 ) 
shows only a small part of the wall being torn away largely intact, 
and (71 = 0.001, 72 = 0 . 006 ) shows the wall brealdng into several 
large pieces. Both of these behaviors demonstrate how the fracture 
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Figure 4: A solid cylinder that experiences ductile fracture when it 
is pulled apart. 



Figure 5: A thin sheet that has been torn apart. 


process can be affected by otherwise unnoticeable amounts of plas- 
tic deformation. The proceedings DVD contains animations that 
further illustrate the behaviors depicted in figure 3 as well as the 
behaviors shown in the other figures. 

Figure 4 shows a solid cylinder tearing as it is pulled and twisted 
apart. Figures 5 and 7 show the ductile fracture that results when 
other objects are ripped apart. 

One way to assess the realism of an animation technique is by 
comparing it with the real world. Figure 6 shows a real clay slab 
that has been struck by a spherical projectile and a simulated slab 
of plastic material that has also be struck by a spherical projectile. 
Although the two images have obvious differences, the holes left 
by the projectiles demonstrate qualitative similarities. 

While modifying the computation of the element stresses to use 
the elastic strain instead of the total strain requires only minor 
changes to an existing code, the change may also have an effect 
on the integration scheme. Our implementation uses an explicit in- 
tegrator that takes adaptive time steps. The step size is determined 
by monitoring the total energy to ensure that the system is not go- 
ing unstable. We compared the size of steps taken when simulating 
a purely elastic material to those taken when simulating a material 
that was identical except that the plasticity code had been enabled. 
During periods when collisions were occurring, both simulations 
took similar-sized integration steps. At other times, however, the 
average step size for the plastic material was approximately twice 
that of the purely elastic one. This result is not surprising because 
plastic deformation absorbs energy implying that it should tend to 
help stabilize the system, but it is only a single test on a single set 
of parameters and further tests would need to be done before any 
more general statement could be made. 

The deformation model we implemented allows a regime of elas- 
tic deformation, followed by a plastic regime, and then possibly 
followed by a second elastic regime. While this model suffices for 
many materials, other materials, such as woven fabrics, may go 
through multiple cycles of elastic and plastic behavior. We have 
also worked only with a linear relationship between elastic strain 
and stress. While a linear model adequately describes many materi- 
als, other materials such as biological tissues demonstrate distinctly 
non-linear elastic behavior. Developing adequate graphical models 
for these types of materials remains an area for future work. 
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Figure 7: A cartoon character being dismembered by a red torture 
device. 
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Abstract 

This paper describes a real-time technique for generating realistic 
and compelling sounds that correspond to the motions of rigid ob- 
jects. By numerically precomputing the shape and frequencies of an 
object’s deformation modes, audio can be synthesized interactively 
directly from the force data generated by a standard rigid-body sim- 
ulation. Using sparse-matrix eigen-decomposition methods, the de- 
formation modes can be computed efficiently even for large meshes. 
This approach allows us to accurately model the sounds generated 
by arbitrarily shaped objects based only on a geometric description 
of the objects and a handful of material parameters. We validate our 
method by comparing results from a simulated set of wind chimes 
to audio measurements taken from a real set. 

CR Categories: 1.3.5 [Computer Graphics]: Computational 

Geometry and Object Modeling — Physically based modeling; 
1.3.7 [Computer Graphics]: Three-Dimensional Graphics and 
Realism — Animation; 1.6.8 [Simulation and Modeling]: Types of 
Simulation — Animation; H.5.5 [Information Interfaces and Presen- 
tation]: Sound and Music Computing — Signal analysis, synthesis, 
and processing 

Keywords: Sound modeling, physically based modeling, sim- 

ulation, surface vibrations, dynamics, animation techniques, finite 
element method, modal synthesis, modal analysis. 


1 Introduction 

One of the central goals for the field of computer graphics is the 
compelling portrayal of realistic synthetic environments. However, 
generating convincing animations of scenes such as that shown in 
figure 1 requires depicting not only the visual aspects of the scene, 
but its audio components as well. While constructing a soundtrack 
by hand often provides a feasible option for animations that are gen- 
erated off line, interactive applications increasingly rely on physi- 
cally based simulation techniques to generate animated motions in 
real-time and these applications require methods for generating the 
corresponding audio in real-time as well. 

One class of simulation method that has found widespread use 
in real-time applications is rigid-body simulations. Because rigid 
bodies are made up of incompliant materials, they experience only 
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Figure 1: A synthetic environment containing a set of simulated 
wind chimes. Both the motion of the chimes and the corresponding 
audio can be computed at interactive speeds. 


small- amplitude deformations during interactions with their envi- 
ronment. Explicitly discarding these small deformations allows 
rigid-body simulators to model a system’s remaining degrees of 
freedom efficiently. However, although visually insignificant, it is 
the vibration of these small- amplitude deformations that generates 
the sounds heard from these objects. 

This paper describes a real-time technique for generating real- 
istic and compelling sounds that correspond to the motions gener- 
ated by rigid-body simulation methods. Precomputing the shape 
and frequencies of an object’s deformation modes allows that ob- 
ject’s vibrational response to contact forces to be efficiently com- 
puted at runtime. The vibrational response is then used directly 
to compute the corresponding audio. Our technique computes an 
object’s deformation modes numerically by performing an eigen- 
decomposition of the system matrices from a finite element model 
of the object. This approach allows us to accurately model the 
sounds generated by arbitrarily shaped objects based on a geomet- 
ric description of the object and a handful of material parameters. 
The diagram in figure 2 provides an overview of this process. 


2 Background 

The technique presented in this paper is closely related to previous 
methods developed by van den Doel, Kry, and Pai. The concept of 
using the vibrational modes of an object for generating sound was 
originally introduced to the graphics community in [van den Doel 
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Figure 2: This diagram illustrates both the preprocessing steps that 
are used to construct the audio/visual model for an object, and the 
processes that subsequently generate sound and motion from this 
description at interactive speeds. 

and Pai, 1996] and [van den Doel and Pai, 1998]. They showed 
that the analytically computed vibrational modes of simply shaped 
objects, such as square plates or cylindrical rods, could be used to 
generate sound maps over the object surfaces. They constructed a 
system that computed realistic contact sounds when an interactive 
user indicated a point on the surfaces of the modeled objects. Be- 
cause the maps were generated from the mode shapes, the system 
correctly captured the variations that arise when objects are struck 
at different locations. More recently, [van den Doel et al., 2001] de- 
scribed a method that uses recorded data to construct sound maps 
over the surface of an object. In addition to allowing a user to inter- 
actively generate sounds by tapping the object surfaces, they used 
contact forces from a real-time rigid-body simulation to excite the 
sampled modes. 

Our work builds on the ideas of these two previous methods by 
extending the range of objects that can be modeled to include ones 
that are neither simple shapes nor available to be measured. The 
analytically computed modes used in [van den Doel and Pai, 1996] 
and [van den Doel and Pai, 1998] are the continuous equivalent 
of the numerically computed discretized modes described in Sec- 
tion 3.1 of this paper. Numerical computation allows determining 
the modes for essentially arbitrary shapes as opposed to a few sim- 
ple shapes, and it makes fewer assumptions about the underling 
differential equations. Our method for driving the sound generation 
from a rigid-body simulation is essentially identical to the method 
used in [van den Doel et al., 2001], but we have found that useful 
results can be generated using “off-the-shelf” rigid-body simulators 
and that a special contact method is not necessarily required unless 
one wishes to produce rubbing or scraping sounds. 

Another related method for generating sound has been described 
in [O’Brien et al., 2001]. Their approach uses a nonlinear finite 
element method to explicitly model the response of an object to 
external forces. Audio is generated by analyzing the computed sur- 
face behavior and then applying a set of filters to the computed 
motion for extracting frequency components that fall within the au- 
dible range. Unlike the method detailed in this paper and the previ- 


ously described methods of van den Doel and his colleagues, the use 
of a nonlinear finite element method allows them to model sounds 
that arise from nonlinear behaviors such as buckling. The main 
limitation of their method is that it requires large amounts of com- 
putation. In contrast, our method can accurately model only sounds 
produced by linear phenomena, but it can compute these sounds in 
real-time. 

In addition to the above physically motivated work on sound gen- 
eration, other prior work in the graphics community has focused 
on sound propagation and heuristic approaches to sound genera- 
tion. Producing synchronized soundtracks for animations was ad- 
dressed in [Takala and Hahn, 1992] and [Hahn et al., 1995]. For 
modeling tearing cloth, [Terzopoulos and Fleischer, 1988] gener- 
ated soundtracks by playing a pre-recorded sound whenever a con- 
nection in a spring mesh failed. Work described in [Savioja et al., 
1997] focused on creating virtual musical performances in virtual 
spaces using physically derived models of musical instruments and 
acoustic ray- tracing for spatialization of the sound sources. Other 
researchers have developed methods for correctly modeling reflec- 
tions and transmissions within the sonic environment [Funkhouser 
et al., 1998; Funkhouser et al., 1999; Min and Funkhouser, 2000; 
Tsingos et al., 2001]. 

The method described in this paper is also related to previous 
work in the graphics community on modeling deformable objects. 
The idea of decoupling an object’s rigid-body behavior from it’s 
elastic deformation was proposed in [Terzopoulos and Fleischer, 
1988] as an efficient method for modeling deformable objects. This 
idea was extended in [Pentland and Williams, 1989] by using modal 
analysis for modeling deformable objects, although instead of ac- 
tually using the object’s vibrational modes they approximated them 
with arbitrary linear and quadratic deformation fields. 

Outside the field of computer graphics, an extensive amount of 
research on sound modeling has been conducted in the digital sound 
and music communities. There the focus has been primarily on ac- 
curately modeling the sounds generated by musical instruments, in- 
cluding the fine subtleties that distinguish high-quality instruments. 
A comprehensive review of the work that has been done in those 
areas can be found in [Cook, 2002]. 

3 Methods 

The mechanical dynamics of a solid physical object can be decom- 
posed into two components: idealized rigid-body motions and elas- 
tic deformations. An object is referred to as being rigid, or incom- 
pliant, if its response to typical interactions includes only negligi- 
ble deformations. For example, when a person taps the side of a 
drinking glass it fiexes slightly but the amplitude of this deforma- 
tion is small enough to be unobservable by sight. However, this 
small deformation may be observable by hearing. In particular, if 
the elastic properties of the glass are such that the small deforma- 
tion induced by the tap results in vibrations at frequencies between 
approximately 20 and 20, 000 Hz, then the small pressure fiuctua- 
tions caused by the oscillating deformation will be heard as sound. 
For further information on the physical process of sound generation 
we refer the reader to [Kinsler et al., 2000]. 

3.1 Modal Analysis 

Our method for modeling the sounds generated by rigid objects 
makes use of a well studied technique known as modal analysis. 
This section presents a brief overview of modal analysis and pro- 
vides the framework for describing our methods. We refer the 
reader to [Cook et al., 1989] for additional information on modal 
analysis. 

A physical system that has been discretized using a finite ele- 
ment, finite differencing, or other similar method can be expressed 
in the following general form: 

K{d) + C{d,d)^Ad(d) = f (1) 
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where d is the vector of node displacements, an overdot indicates 
a derivative with respect to time, K and C are nonlinear functions 
that respectively determine the internal forces due to node displace- 
ments and node velocities, A1 maps node accelerations to node 
momenta, and f represents any other (e.g. external) forces. Typi- 
cally, the forces determined by K are internal elastic forces and C 
determines damping forces. 

In general, equation (1) is nonlinear, however if we assume that 
the displacements are small then we may linearize about the sys- 
tem’s rest configuration giving: 

Kd + Cd + Md = f (2) 

where K, C, and M are respectively known as the system’s stiff- 
ness, damping, and mass matrices. For the physical systems corre- 
sponding to solid objects, all three matrices are real and symmet- 
ric. Both K and C are positive semi-definite, and M is positive 
definite. Linearizing in this fashion is consistent with our goal of 
modeling the small- amplitude, high-frequency vibrations in solid 
objects that produce sound. Unfortunately, the linearized system 
cannot model the rotational components of rigid-body motion. We 
will put this issue aside for now, but later we will return to it and 
show how the rigid-body modes can be decoupled from all other 
modes. 

Once we have the linearized system, the next step in the modal 
analysis is to perform a series of manipulations that will diagonal- 
ize equation (2). To facilitate this process, we will first assume 
that C = aiK a 2 M for some ai and 0 ^ 2 . Expressing the 
damping matrix as a linear combination of the stiffness and mass 
matrices is known as Raleigh damping. Although this assumption 
simplifies diagonalization while still producing good results, it is 
not strictly necessary. A more general assumption, known as pro- 
portional damping, that expresses the damping matrix as a linear 
combination of powers of the stiffness and mass matrices would 
also be diagonalized by the process described below but the equa- 
tions would be more cumbersome. Additionally, even if for some 
reason C must be arbitrary, then other, slightly more complicated, 
methods are available for decoupling equation (2) [Anderson et al., 
1999; Bai et al., 2000]. 

Replacing C with aiK a 2 M gives: 

K{d+aid) + M{a2d + d) = / . (3) 

Since M is symmetric and positive definite, it may be decomposed 
using a Cholesky factorization so that M = LL^ . If we introduce 
another variable, y — d, and then rewrite equation (3) in terms 
of y after pre-multiplying hy L~^ we then have: 

L-^KL-^{y + a^y) + {a2y + y)=L-^f . ( 4 ) 

The real and symmetric matrix L~^KL~^ can be decomposed 
into L-^KL-^ = VKV^ where V is the orthogonal ma- 
trix whose columns are the eigenvectors of L~^KL~^ and A is 
the diagonal matrix of eigenvalues. Introducing another variable, 

= V^y, and pre-multiplying by transforms equation (4) 
into: 

A{z + aiz) + {a 2 Z + z) = V'^ L-^ f (5) 

which can be rearranged to give: 

Az + (aiA + a 2 l)z z = g (6) 

where g = V'^L-^f. 

At this point the original linear system of equation (3) has been 
diagonalized into a set of decoupled oscillators. The z’th row of 
equation (6) is the scalar second-order differential equation: 

XiZi + {aiXi + a2)zi + Zi — gi ( 7 ) 


where Xi is the z’th entry of the diagonal matrix A. Equation (7) 
may be solved by numerical integration or it may be solved more 
efficiently using the analytic solution: 

tu}~y I tuj . /o\ 

Zi = cie ^ + C2e ^ (8) 

where ci and C 2 are arbitrary (complex) constants, and Ui is the 
complex frequency given by 

—{aiXi + a2) i (o^iA^ + 0 ^ 2 )^ — 4Ai 

^ ^ 

The absolute value of the imaginary part of cui is the frequency 
(in radians/second, not Hertz) of the mode, and the real part is the 
mode’s decay rate. 

The decoupled system of equation (6) is not an approximation 
of the original linear system in equation (3), it is exactly the same 
as the original linear system. Of course the linear system was an 
approximation of the original nonlinear one, but any problem that 
could be solved using equation (3) could also be solved with equa- 
tion (6). 

The columns of L~^V are the vibrational modes of the object 
being modeled. (See figure 3.) Each mode has the property that 
a displacement or velocity over the object that is a scalar multi- 
ple of the mode will produce an acceleration that is also a scalar 
multiple of the mode. This property means that the modes do not 
interact with each other, which is why decoupling the system into a 
set of independent oscillators was possible. The eigenvalue for each 
mode is the ratio of the mode’s elastic stiffness to the mode’s mass, 
and it is the square of the mode’s natural frequency (in radians per 
second). In general the eigenvalues will be nonzero, but for each 
free body in the system there will be six zero eigenvalues that cor- 
respond to the body’s six rigid-body modes. The rigid-body eigen- 
values are zero because a rigid-body displacement will not generate 
any elastic forces. 

3.2 Rigid Body Simuiation 

As discussed previously, the rigid-body modes for an object do not 
interact with the object’s deformation modes provided the amount 
of elastic deformation experienced by the object is small. ^ Addi- 
tionally, small- amplitude elastic deformations will not significantly 
effect the rigid-body collisions between objects. These observa- 
tions allow us to model the rigid-body behavior of the objects in al- 
most the same way as if we were not interested in generating audio. 
The only change that must be made to the rigid-body simulation is 
that information about contact forces must be gathered and exported 
to another process that will generate the audio. Of course, hearing 
the results of the rigid-body simulation, in addition to seeing them, 
may reveal previously unnoticed inadequacies of the simulator, but 
we have not found this to be a problem with the simulation engines 
we have worked with. 

We have implemented our system using two existing rigid-body 
simulation engines that were not originally designed for generating 
audio. Our choice of engines was motivated by what systems were 
readily available and how well they were able to model the scenar- 
ios we wished to test. The first is a commercial software package. 
Vortex, sold by Critical Mass Labs. The second system we are us- 
ing had been previously written by Okan Arikan, a graduate student 

^Actually, the requirement was that all displacements be small, includ- 
ing displacements corresponding to the rigid-body modes. The translation 
modes are inherently linear so they cannot interact with the elastic modes 
regardless of their magnitude, but for a rapidly rotating body there will be 
some coupling between the rotation modes and the elastic ones. Unless the 
object is rotating very rapidly or experiencing large angular accelerations, 
the coupling between rotation and elastic modes with frequencies in the au- 
dible range will be negligible, so we ignore this interaction. 
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not involved in this project. No special changes were made to either 
package other then instrumenting them to allow reporting collision 
forces. 

3.3 Deformation Model 

Once the task of modeling the rigid-body modes has been dele- 
gated to a rigid-body simulator, the remaining elastic deformation 
modes can be used for generating audio. Because we are interested 
in modeling sounds from incompliant objects, we can use the modal 
decomposition methods described in Section 3.1 to compute their 
behavior efficiently. However before we can perform a modal de- 
composition, we must first select a deformable modeling method 
that can be used to generate the K, C, and M matrices. 

The method we are using for modeling deformable behavior 
is the tetrahedral finite element method described by [O’Brien 
and Hodgins, 1999] for modeling fracture propagation, and subse- 
quently used in [O’Brien et al., 2001] for modeling nonlinear audio 
generation. As discussed by O’Brien, Cook, and Essl, a variety of 
methods could be used, including spring/mass systems or finite dif- 
ferences methods. We selected this finite element method because 
their previous results show that it is accurate enough for generating 
compelling audio. 

Computing the global stiffness and mass matrices proceeds by 
first computing individual 12x12 stiffness and mass matrices for 
each element and then assembling the results to form the global 
matrices. From [O’Brien and Hodgins, 1999] the nonlinear node 
forces are given by: 


/[i]a = 


kCTkl 


( 10 ) 


i=i 


where f^-^a is the a’th component of the force exerted on the z’th 
node of the element, vol is the volume of the element, p are the 
node positions, (3 is the element basis matrix, and a is the stress 
tensor within the element. Details for computing (3 and a appear 
in [O’Brien and Hodgins, 1999]. 

The element stiffness matrix, fc, is computed by taking the par- 
tials of f and evaluating them at zero displacement: 


\ij]ah 


df[i\a 


9p[j]b 


( 11 ) 


P=Pr. 


— 2 ^ f^ikf^jkSgb) ( 12 ) 


where 6 is the Kronecker delta, and A and p are the material’s Lame 
constants.^ This is the exactly the same matrix that would have re- 
sulted if Cauchy’s infinitesimal strain had been used in place of 
Green’s strain, however with Cauchy’s strain the partials would be 
constant with respect to node position so that it would not matter 
where they were evaluated. 

The element mass matrix, m is computed by taking the second 
partials of the kinetic energy within the element with respect to the 
node velocities, which turns out to be constant with respect to node 
position and velocity: 


'fT^[ij]ab 


^P[i]a^P[j]b 

^ii+Sij)Sab 


(13) 

(14) 


where k is the kinetic energy within the element, an overdot repre- 
sents a derivative with respect to time (i.e. p are node velocities), 
and p is the material’s density. 

The global stiffness and mass matrices, K and iVf, are built by 
assembling the element matrices. Assuming that we are working 
with three-dimensional objects, each of the global matrices will be 
3N X 3N where N is the number of nodes in the finite element 
mesh. Each entry in each of the 12x12 element matrices is accu- 
mulated into the corresponding entry of the global matrix. 

Since each node in a tetrahedral mesh will share an element with 
only a small number of the other nodes, the global matrices will be 
very sparse. This sparseness means that an eigen decomposition of 
K can be performed efficiently using sparse matrix algorithms. Un- 
fortunately, the Cholesky decomposition tends to generate a dense 
L matrix even when M is originally sparse, and as a result com- 
puting L~^ may be costly and L~^KL~^ will be densified. 

Dense matrix algorithms can be used for systems up to approx- 
imately 1000 nodes, but beyond that we suggest using an alternate 
mass matrix that does not generate a dense Cholesky decomposi- 
tion. The alternate mass matrix, known as a lumped mass matrix, 
simply shifts the sum of each row onto the diagonal: 


lumped pvol 

^[ij]ab - 


3ij Sab • 


(15) 


Because the element mass matrices are diagonal, the global mass 
matrix will be as well, and its Cholesky decomposition will also 
be diagonal: L will be a diagonal matrix whose entries are simply 
the square root of the entries of the lumped M. For small systems 
generated by coarse meshes, the errors introduced by mass lump- 
ing may be significant. However, as the mesh gets finer the errors 
introduced by lumping quickly become insignificant [Cook et al., 
1989]. Luckily, the large systems corresponding to fine meshes are 
precisely the ones that require the sparse solvers facilitated by mass 
lumping. Our implementation includes both dense and sparse de- 
composition routines and we use whichever is appropriate to the 
size of a particular system. For dense decompositions, we use the 
routines from LAPACK [Anderson et al., 1999], and for sparse de- 
compositions we use the TRLan package [Wu and Simon, 1999]. 
The method used for each of our examples, along with computation 
times and the number of nodes, is listed in table 1. 

The use of Raleigh damping was another simplification that we 
made to facilitate decoupling equation (2). In [O’Brien and Hod- 
gins, 1999] they used a nonlinear stiffness-proportional damping 
term based on the strain rate with parameters 0 and V’. Raleigh 
damping is equivalent to a linearization of this damping term with 
the additional constraint that ^ and the Raleigh parameter ai 
should be set to this ratio to generate equivalent results. O’Brien 
and Hodgins did not discuss a mass proportional damping term, but 
setting 0(2 to a non-zero value would be equivalent to including a 
{—a 2 dimi) damping force on each node. 

Even with sparse matrix methods, computing a system decompo- 
sition still requires a significant amount of time, so it is worth noting 
that certain changes may be made without recomputing the decom- 
position. The damping parameters, oi and 02 , have no effect on 
the decomposition, so the only work involved when changing them 
is re-evaluating equation (9). Changing the material’s density does 
not change the mode shapes, it only scales the eigenvalues by the 
inverse of the scale factor applied to the density. Similarly, scal- 
ing the Lame constants both by the same scale factor (i.e. so that 
the ratio between A and p is preserved) only scales the eigenvalues 
by the same ratio. Changing the ratio between the Lame constants, 
changing the shape of the object, or modifying the mesh all require 
recomputing the decomposition. 


^Unfortunately, the symbol A is conventionally used both to indicate one 
of the system eigenvalues and the first Lame constant. In this paper it should 
be clear from context (and the presence or absence of a subscript) what the 
symbol is referring to. 


3.4 Sound Generation 

Once all of the computational machinery described above is avail- 
able, the actual process of computing audio matching the motion 
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Figure 3: The top shows a multi-exposure image from an anima- 
tion of a bowl falling onto a hard surface with the path of the bowl’s 
center traced by a yellow curve. Only the bowl is sounding. The 
two bottom rows show a side and top view of the bowl along with 
three of the bowl’s first vibrational modes. (The modes selected for 
the illustration are the first three non-rigid ones with distinct eigen- 
values that are excited by a transverse impulse to the bowl’s rim.) 


from a rigid-body simulation is both straightforward to implement 
and computationally efficient: 

1. A rigid-body simulation is set up for the desired scenario. 


modes with |Im(a;i)| in the range 3.18 ... 3,180 Rad/s are used, 
and then simplifying yields 

sin(i|Im(ci;i)|) . ( 18 ) 

j I 

Evaluating equation (18) for every audio sample is inefficient. 
By noting that the value of the oscillator at one 

audio sample can be computed from the previous value using only a 
single complex multiply. Additionally, as a mode is excited at sub- 
sequent times by different contact forces, the additional excitations 
can be modeled by simply adding the new value to the oscillator’s 
current value. Because the cost of modeling additional impulses is 
essentially zero, the forces from the rigid-body simulation may be 
convolved with a Gaussian kernel to model the effect of soft colli- 
sions, or with a noise function to model small-scale roughness that 
is below the resolution of the rigid-body simulator [van den Doel 
et al., 2001]. Our results were generated using the former. 

A method for modeling the coupling between vibrations in an 
object and vibrations in the surrounding air is described in [O’Brien 
et al., 2001]. Unfortunately, their method is too slow for real-time 
use. We compute an approximate coupling coefficient for each 
mode by summing the amount of normal displacement generated by 
that mode over the surface of the object multiplied by the mode’s 
frequency. The coupling coefficient for each mode multiplies the 
result computed by that mode’s oscillator and the sum of the scaled 
oscillators is the final sound generated by the system. A result 
of this simplification is all objects are treated as omni-directional 
sources. 


4 Results and Discussion 


2. For each object in the simulation, the system matrices are as- 
sembled and decomposed into their vibrational modes (i.e. the 
columns of L~^V). 

3. For each object in the simulation, only the columns of L~^V 
corresponding to |Im(ct;i)| in the range 3.18 ... 3,180 Rad/s 
(20... 20, 000 Hz) are retained, the rest are discarded (or if the 
sparse method is used, never computed). 

4. As the rigid body simulation runs, collision forces are pro- 
jected onto the retained modes. The response of each mode is 
modeled using equation (8). 

5. Each mode response is scaled according to how it moves the 
objects surface and the scaled responses are then summed to- 
gether. 

6. Einally, the result is output to the computer’s audio device. 


In practice, not all of the modes in the audible range need to 
be retained. As discussed in [van den Doel et al., 2001], high- 
quality results can easily be obtained using only the first 800 or 
fewer modes. 

A mode’s response to a projected impulse is given by equa- 
tion (8) with 


2Atgi 

Cl = ^ 3 

2Atgi 

C2 = ^ + 

^z 


(16) 

(17) 


where At is the interval over which the projected force is applied, 
and t is time relative to when the impulse was applied. Substitut- 
ing these values of ci and C2 into equation (8), recalling that only 


We have built a system that implements the methods described 
above and used it to generate a number of demonstrative examples. 
Table 1 lists the parameters that were used in each of the examples, 
and the video tape accompanying this paper contains animations 
that exhibit the sounds and motions produced. 

To test how well the computed results match real objects, we 
generated the wind chimes shown in figure 1. These chimes were 
modeled based on measurements from a real set of chimes. Each 
tube is a hollow cylinder 1.25 cm in radius with a nominal wall 
thickness of 1 mm. The measured lengths of the chimes are listed 
in table 2. We computed the modal decomposition for each chime 
using reference parameters for aluminum. The resulting base fre- 
quencies matched measured ones to within 2% error. However, the 
real chimes were slightly out of tune, so we tuned the simulated set 
by adjusting the tube lengths so that they were within ±1 Hz of the 
correct (D scale) tuning. 

Eigure 3 shows a bowl model that was used for two of the ex- 
amples. The modal decomposition of the bowl was computed once 
with material parameters for aluminum and again with material pa- 
rameters for wood (oak). Two animations were created, both with 
the same rigid-body motion but with the two sound tracks gener- 
ated from the two different modal decompositions. The resulting 
audio (refer to video tape) captures the general characteristics of 
both materials as well as details such as the sound produced as the 
bowl rolls on its edge. Eigure 3 also illustrates the mode- shapes 
for three of the bowl’s vibrational modes by showing the results of 
applying the mode as a displacement over the bowl’s original shape. 

An example generated using a more complex model consists of 
bunny figurines falling through a chute. (See figure 4.) Both the 
bunny and the shelves in the chute generate sounds when struck. 
The shelves are made of plastic, metal, and wood. The bunny is 
ceramic. The tetrahedral bunny model was generated by meshing 
the region between the surface of the Stanford Bunny model and an 
interior offset surface to create a hollow figure with finite thickness 
walls, as shown on the right side of figure 4. The right side of 
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Example 

Eigure 

A (Pa) 

M (Ps-) 

ai 

a2 

p (Kg/m3) 

Base Ereq. (Hz) 

Decay 

Num. Nodes 

Method 

Precompute 

Chime(D3) 

1 

4.98 X 10^° 

2.57 X 10^® 

1 X lO-"^ 

0 

2700 

587.4 

0.6 

18796 

Sparse 

2h 24min 

Bowl #1 

3 

4.98 X 10® 

2.57 X 10® 

1 X lO-"^ 

30 

2700 

551.3 

15.6 

387 

Dense 

4min 12sec 

Bowl #2 

- 

5.00 X 10® 

1.00 X 10® 

8 X 10“® 

50 

750 

216.7 

22.4 

387 

Dense 

4min 12sec 

Bunny (Ceramic) 

4 

3.99 X 10® 

2.05 X 10® 

1 X 10“® 

10 

2700 

855.9 

19.5 

37114 

Sparse 

4h 40min 

Plastic Shelf 

4 

2.49 X 10^° 

1.28 X 10^® 

1 X 10“® 

50 

2700 

488.9 

29.7 

361 

Sparse 

30sec 

Aluminum Shelf 

4 

4.98 X 10^® 

2.56 X 10^® 

1 X lO-"^ 

0 

2700 

691.5 

0.9 

361 

Sparse 

30sec 

Wood Shelf 

4 

5.00 X 10® 

1.00 X 10® 

8 X 10“® 

50 

750 

154.6 

28.8 

361 

Sparse 

30sec 

Bunny (Metal) 

- 

4.99 X 10^® 

2.56 X 10^® 

1 X lO-"^ 

0 

2700 

855.9 

19.5 

37114 

Sparse 

4h 40min 

Blocks 

5 

5.00 X 10® 

1.00 X 10® 

8 X 10“® 

50 

550 

1596.2 

428.1 

1160 

Dense 

5h 28min 

Boxes 

5 

5.00 X 10® 

1.00 X 10® 

8 X 10“® 

50 

550 

159.1 

49.0 

1160 

Dense 

5h 28min 

The End (T) 

6 

1.49 X 10® 

7.70 X 10® 

2 X lO-"^ 

30 

2700 

247.7 

15.2 

71 

Dense 

42sec 


Table 1: This table lists parameters that were used for each example object, the resulting frequency and decay for the object’s primary mode, 
the number of nodes in the tetrahedral mesh, the method used for the modal decomposition, and the amount of time required to compute the 
decomposition. Once the model decomposition has been computed, all of the above examples can generate audio in real-time. For “Chimes” 
and “The End,” the information listed is for the D3 tube and the letter T. 


Note 

Ideal 

Ereq. 

Measured 
Length Ereq. 

Computed 

Ereq. 

Adjusted 

Length Ereq. 

D3 

587.33 

.505 

585.8 

589.17 

.5061 

587.40 

E3 

659.26 

.475 

656.0 

665.03 

.4770 

659.27 

G3 

783.99 

.435 

781.8 

787.01 

.4366 

784.06 

A4 

880.00 

.410 

877.5 

884.70 

.4115 

879.36 

B4 

987.77 

.388 

982.5 

984.75 

.3878 

987.32 

D4 

1174.66 

.353 

1167.0 

1186.88 

.3548 

1174.67 


Table 2: The notes and ideal frequencies listed indicate the values 
specified by the manufacturer of the real wind chimes. The mea- 
sured values were taken from the real wind chimes. The computed 
frequencies are what our model produced using the parameters from 
table 1 and the measured lengths. The adjusted values indicate the 
length and resulting frequency of the simulated chimes after tuning. 
Lengths are in meters and frequencies in Hertz. 

figure 4 also shows the results of projecting a pair of impulses onto 
the retained modes of the bunny model. 

The blocks and boxes shown in figure 5 illustrate how scale can 
effect the resulting audio. Both the boxes and blocks are geometri- 
cally similar: hollow cubes with a wall thickness of 5% their width. 
However, the boxes are 10 x the size of the blocks. While the dif- 
ferent scales are subtly revealed by the rigid-body motions (by the 
rate of acceleration with respect to the object sizes), the sounds 
produced by the two sets of objects are distinctly different, and the 
difference provides a clear cue as to the size of the objects. 

As we discussed previously, similarities exist between the ap- 
proach we have presented here and that presented in [van den Doel 
et al., 2001]. The main difference between the two methods is that 
we synthesize audio from only geometry and material properties 
whereas their system makes use of extensive measurements of a 
given object’s response to impacts. Each of these methods presents 
distinct advantages: by relying on recorded data their method may 
easily match a given object, but our method is applicable when no 
real object or no robotic measuring devices are available. One di- 
rection that might be worth pursuing would be using their measured 
data for a given object to infer material parameters that could then 
be applied to the geometry of a different object. This approach 
might allow audio models for an entire set of cooking pots, for ex- 
ample, to be generated from measurements of a single pot in the 
set. It might also allow us to determine the sound made by a novel 
bell design, based on data from bells of similar materials, before we 
actually make the bell. Based on the good correspondence between 
our synthetic chimes and the physical set, we are optimistic about 
this direction of future work. 



Eigure 4: The left side of this figure shows an image from ani- 
mation of several bunnies falling through a chute. Both the bunnies 
and the shelves are sounding. The images on the right show in order 
from top to bottom: the exterior of the bunny model, a cut-away re- 
vealing the wall thickness and hollow interior, modal response to an 
impulse on the bunny’s nose, and the modal response to an impulse 
on the bunny’s back. The impulse responses are greatly exaggerated 
for illustration. 


Although the resolution of the mesh can affect the resulting au- 
dio, we have found that even very coarse meshes may be used for 
generating acceptable results. The meshes used for each of the let- 
ters shown in figure 6 are very coarse, yet the resulting audio is still 
acceptable. We have found that low mesh resolution tends to shift 
frequencies higher and may add a “hollow” quality to the sound. 
The frequency shifting may be partially compensated for by simply 
modifying the material parameters (e.g. raising the density) to com- 
pensate, so it will only be a problem if one is attempting to match a 
particular object (as we were for the wind chimes). 

Although the modal decompositions may require up to a few 
hours of computation, this work needs only to be done once for 
a given object and audio can then be generated interactively. By 
precomputing the modal decomposition and storing it with an ob- 
ject, the approach we have presented could easily be applied to 
interactive applications such as video games that already employ 
rigid-body simulation methods. Additionally, because our method 
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Figure 5: The top images show a stack of 5 cm blocks being 
knocked over. The images on the bottom show a stack of 50 cm 
boxes being knocked over. Only the blocks and boxes are sound- 
ing. Other than the 10 x scale, both models are identical. The plots 
below each sequence show the frequency content of the resulting 
audio, indicating a significant difference in the sounds. (The hor- 
izontal axis ranges from 0 to 5000 Hz, the vertical axes are auto- 
scaled independently.) 


requires only a geometric model and a handful of material param- 
eters, the extra effort required to generated the audio model of a 
given object is minimal. 
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Figure 6: Select frames from an animation of the words “The End” 
falling onto a hard surface. Both the letters and the surface are 
sounding. 
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